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ABSTRACT 

A two-dimensional hydrodynamics code for Type la supernovae (SNIa) simulations is 
presented. The code includes a fifth-order shock-capturing scheme WENO, detailed 
nuclear reaction network, flame-capturing scheme and sub-grid turbulence. For post¬ 
processing we have developed a tracer particle scheme to record the thermodynamical 
history of the fluid elements. We also present a one-dimensional radiative transfer code 
for computing observational signals. The code solves the Lagrangian hydrodynamics 
and moment-integrated radiative transfer equations. A local ionization scheme and 
composition dependent opacity are included. Various verification tests are presented, 
including standard benchmark tests in one and two dimensions. SNIa models using 
the pure turbulent deflagration model and the delayed-detonation transition model 
are studied. The results are consistent with those in the literature. We compute the 
detailed chemical evolution using the tracer particles’ histories, and we construct cor¬ 
responding bolometric light curves from the hydrodynamics results. We also use a 
Graphics Processing Unit (GPU) to speed up the computation of some highly repet¬ 
itive subroutines. We achieve an acceleration of 50 times for some subroutines and a 
factor of 6 in the global run time. 


1 INTRODUCTION 

1.1 Chandrasekhar Mass Explosion Model 


A Type la supernova (SNIa) is the explosion of a carbon- 
oxygen white dwarf (WD) due to thermonuclear runaway of 
carbon burning. It is believed to be a s tandard candle due to 
its o bserved explosion homogeneity dBranch fc Tammann . 
1992) and the luminosity-width relation (IPhillips et al. , 
19871 ). Also, in numerical modeling, the standard Chan¬ 


drasekhar mass WD is regarded as the progenitor of an 
SNIa. These properties lead to wide applications of SNIa 
as a cosmological ruler for deter mining the Hubble parame¬ 
ters (|Leibundgut&Pintd, 1992|), and the d i scove ry of dark 
energy ( Riess et al. . 1 19981 : IPerlmutter et all 1 19991 ) . 

In the last few decades, several explosion mechanisms 
have been proposed, starting from the pure det onation 
model of Chandrasekhar mass WD proposed in lArnettl 
(119691) . Despite its ability in unbinding the whole star, the 
over-production of iron-peaked elements and the absence of 
intermedi ate mass e lements (IME) make this scheme im¬ 
plausible (|Nomoto fc Sugimotq . Il977l) . Later, the sub-sonic 
pure deflagration model is proposed, which aims at provid¬ 
ing suf ficient time for electron capture in the laminar flame 
stage (iNomoto et al.l . 119761 ) and production of IME. The 
weakness of pure detonation model is resolved but the de¬ 
flagration wave is too slow to unbind the WD. The flame 
is quenched before consuming the whole WD l _wliich leaves 
a lar ge amount of unburnt material (fuel) llNomoto et all 
119761) . In view of this dilemma, several flame acceleration 
schemes or transition schemes are proposed. Popular mod¬ 
els include the delayed-detonation transition model (DDT), 


the pure turbulent deflagration (PTD) model and the grav¬ 
itationally confined detonation (GCD) model. 


The D DT mod el is first sugg ested by Khokhlov (see 
for example IKhokhlov! (Il989l . fl991a U). It is believed that the 
eddies around the flame front may provide the required envi¬ 
ronment, which seed s th e d eton ation spot by the Zel’dovich 
gradient mechanism (lKhokhlovlll991bl) . It has a counterpart 
in the shocktube experiment. The model has been found 
satisfactory because of the sufficient production of IME, 
abse nce of remnant and consumption of fuel a round the 
core (lKhokhlovlll989| : lGamezo et all|2004l . [20051) . However, 
whether the first detonation spot can be seeded is still an 
open question. It is believed that the velo city fluctuations 
are adequate for triggering the detonation (IKhokhlov et al .1 . 
1995). On the other ha nd, several turbu le nt flame simu¬ 
lation s (see for example Nieme yeil (1 19991 ) . Ilmshenik et alJ 
Jl999l) . and lLisewski et al.l 1 200(1) ) have shown that the de¬ 
veloped velocity fluctuations are far below the required level. 
Also, the detonation front is found to be unsta ble when en¬ 
countering obstacles dMaier fc Niemeveil . l2006l) . This makes 
the robustness of this mechanism in doubt. 


The PTD model assumes that the the fluid is 
highly turbu lent in view of the high Reynolds number 
dNiemever fc Hillebrandtl . Il995bl) . Local velocity fluctua¬ 
tions are important for the flame evolution because the flame 
is constantly stretched and perturbed by t he flui d motion , 
which enlarges the local burning surface ((Khokhlovl . Il994l ) . 
This increases the fuel consumption rate and hence makes 
the flame propagate faster , compar ed with a laminar flame 
under the same condition dr imme 3. Il992l ). Two models are 
commonly considered in turbulent flame modeling. The first 
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model is proposed in (iDamkoehleil . 11939h . It is observed in 
Bunsen flame experiments that the turbulent flame propa¬ 
gation speed Uturb varies wi th local turb ulenc e fluctuations 
v: Uturb = v dNiemever fc Hillebrandtl . Il995al l. The secon d 
mode l is obtained from a theoretical analysis (IPocheaul . 

11994 1. It is found that vtmb = viamVl + C(v/v i am ) 2 
with uiam being the laminar flame speed. The applica¬ 
tions of turbulent flame i n the SNIa explosion s cenario 
are first proposed in lNiemever fc Hillebrandtl dl995a[} . based 
on th e sub-grid turbulence model reported in Clementl 
Jl993if . The model has been found succe ssful in provid¬ 


ing a h ealthy explosion in two-dimensional (Reinecke_e£_al 


2002a) and three-dimensional simul a tions (Rcincckc et al.. 
2002b; iRoepke fc Hillebrandtl . 120051 : iRoepkel .' 20051 ). How¬ 


ever, several drawbacks have been found in detailed three di¬ 
mensional studies. They include the presence of low-velocity 
fuels (ca rbon and oxyg en) in the ejecta and underproduction 
of IME (IRoepke et all 120071) . 

The GCD model, also know n as the Defla g ration 
Failed Detonati on (D FD) model llPlewa fc Kasenl . 120071 : 
iKasen fc Plewai . 120071 ). is similar to the DDT model, which 
starts from the deflagration phase. But the burning is weak 
so that the star remains bound. The hot and burnt ma¬ 
terial causes the WD to expand, floating the hot flame 
by buoyancy to the surface. After then, a large amount of 
fuel remains unburnt. The flame propagat es throughout the 
WD by a surface flow (Ijorden et all 120081 ). The flow finally 
merges at one po int, usually at the point opposite to the 
breakout location (iMeakin et al.l . l2009l ). The converging flow 
heats up the cold and low-density fuel, pushing it deep inside 
the star. Once the squeezed fuel becomes the hot spot for 
detonation as its density exceeds the threshold, a detona¬ 
tion front forms, which burn s the remainin g stellar material 
and unbinds the whole WD (Ijorden et alll2012l) . The qual¬ 
itative difference between DDT model and GCD model is 
that, in the former, the detonation starts from the inside of 
the WD and sweeps outward, while it is the opposite in the 
latter. 

Recent developments of hydrodynamics mo dels have 
consid ered extensively all three models. In iLong et al] 
( 201^), three-dimensio nal PTD models are studied. In 
ISeitenzahl et al.l (1 20121 ). the DDT model is sh own to be able 
to exp lain c ertain SNIa ob se rvatio n data. In iBlondin et al] 
(1201 If) and IBlondin et al.l (120121 ). DDT models in one- 
dimension and two-dimension are studied with their syn¬ 
thetic spectra and light curves. Such fine details provide 
a me an to constrain the explosion model dDessart et all 
120131) . 


1.2 Radiative Transfer 

The modeling of post-explosion light curve and spectrum are 
crucial for discriminating the validity of an explosion mecha¬ 
nism. The governing physics of SNIa light curves are believed 
to be the decay of synthesized 56 Ni into 56 Co in early time , 
and decay of 56 Co into 56 Fe in late time dColgate fc Whitel . 
1969). Simplifi ed models, in c ludin g those with diffusion ap¬ 
proximation dColgate et aid . [l980), or analytic mo dels as¬ 
sumi ng a fireball undergoing homologous expansion dArnetd . 
1 19821) . can already capture primary features of SNIa light 
curves and match a number of observed SNIa. 

Despite that the light curves are well fitted by these 


models, in order to understand the variety of SNIa ob¬ 
servations, as well as their relations with explosion mod¬ 
els, detailed radiative transfer for both b olometric and 
multiple wavebands are important dZhang fc Sutherland! 
Il994l ). The actual problem of radiative transfer can be 
highly nontrivial due to its interactions with hydrodyan- 
mics, the presence of millions to billions of atomic tran¬ 
sition lines, ionization and excitation of atoms, and the 
differential-integral structure of the radiative transfer equa¬ 
tions. These features have been studied in details in the 
last few decades. Physics a nd numerical factors, such as 
the effects of line blanketing dHillierl . [l^OljHillier^&MilleiJ, 
Il99gl) . consistent bound ary conditions dSauer et aid . I2006I ). 
treatment i n opac i ties ( Hoeflich__eT_alJ, 19931 ). acceleration 


techniques dHillierl . 1 199C : iLuc i|200l|) and so on have been 
studied extensively. A number of numerical codes for solv¬ 
ing this problem have been under co nstant development in 
the past few decades (see STELLA dBlinnikov et all 120081 : 
[Blinnikov&J; orokinal. 200d: Blinnikov et al.f 20061 ). ARTIS 
( LucJ" 20051: Sim , 20071 : Krorner fc Siml, 20091). SEP ON A 
( Kasenj^ 2006|]and P HOENIX (lHau schildt fc BarorJ . 120101 : 


ISeelmann et al.l . l2010l : lHauschildt fc Baron . 20 111 ) for the i 

strument papers and realizations). 

Due to the stringent demand in computational re¬ 
source for calculating radiative transfer, in previous stud¬ 
ies, the explosion phase and its obser vational conse¬ 
quen ce s are usually mode le d sep a rately dNomoto et af ~ 


19861: Zhang fc Sutherland! 1 1994 ; iHoeflich et al. . 1995 : 
Nugent et aid . 19971 ). where the radiative transfer part makes 


use of the hydrodyanmics results from some benchmark 
runs. Recent development of computational power allows 
the spectral synthesis to be coupled in the hydrodynam¬ 
ics to form a pipeline following the explosion and homol¬ 
ogous expansion phases. Different explosion mechanisms 
have b een studied with fine details , such as the PTD 
model dFink et ahl. 120141; ILong et all. 20141). DDT model 
dBlondin et all 201 it ISeitenzahl et al. ~ |20JJ2 |) and GCD 
model ( Plewa fc Kasenlj2007l : IKasen fc Plewai . 120071) . 


1.3 Flame Capturing Scheme 

The flame capturing scheme is a technique of describing dis¬ 
continuities with details in sub-grid scales. Such discontinu¬ 
ities can be in forms of fluid type, composition or characters. 
In SNIa simulations, it is essential because the typical width 
of deflagration waves is in the order of centimeters or even 
smaller. Also, their propagation speed is much slower than 
the fluid speed of sound. Therefore, within one time step 
under the Courant-Fredrich-Lewy condition, only a partial 
amount of a fluid element in a grid is burnt and so it is im¬ 
portant to determine the ash-fuel interface in partially burnt 
grids. 

To account for such discontinuities, three types of inter¬ 
face tracking algorith ms are commonly used. The first one is 
the level-set method (lOsher fc Sethianl . ll988l ). This method 
suggests that the discontinuity can be described by the zero 
contour of a scalar function, which is interpreted as the dis¬ 
tance function. This method has advantages in its easy im¬ 
plementation and direct coupling with hydrodynamics. Also, 
the scheme can be generalized directly to arbitrary dimen¬ 
sions. The topological changes are handled naturally wit h- 
out the need of extra spotting mechanism dSethianl . 1200ll ). 
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However, it is al so known to have poor volume conservation 
jRider fc Kothd . Il995l l. The distance function requires reini¬ 
tialization at each iterat ion step in o rde r to preserve the dis¬ 
tance function validity (ISussma n et al.. 1994). Th is scheme 
in SNIa is first propos ed bvlReinecke et al.l 1 1999|L and ap¬ 
plications are found in lReinecke et alT i 19991 l2002al lbh. 

Another alg orithm is the point-set m ethod 
jGlimm_et_alJ, _ 198ll : Idimm fe McBrvanl 1 19851 : 


IGlimm et all Il98ab The flame front is described 
line in a two-dimensional simulation or a surface in a 
three-dimensional simulation, by a linked set of massless 
particles which are advected along the streamline of fluid 
flow, and they only reveal the position of the interface 
without affecting the hydrodynamics properties of the fluid. 
This method receives wide applications in other fields , 
such as the study of bubbles dYoungren fe Acrivosl . Il976l) . 
It has a huge advantage that the interface properties are 
exactly modeled, such that their geome tric properties can 
be ob tained in a straightforward manner dTryggvason et all 
l200ll l. However, it has two major shortcomings. First, it 
requires special attention in t opological changes, namely 
surface splitting and merging (iGlimm et al l. 1 19881 ). Also, 
the extension of the point-set method to three-dimensional 
simulations is non-trivial, becaus e the surface i s no longer 
represented by a line segment JGlimm et all Il999t l . For 
example of application, see IZhang ( 2009l l for a laminar 
deflagration model. It is shown that, with sufficient fine 
details, the laminar flame can bring a successful explosion 
without the need of any flame acceleration scheme. 

The third mechanism which is commonly foun d in the 
litera ture is the volume of fluid method dHirt fe Nicholsl . 
Il98lh . This method is similar to the level-set method, which 
introduces an extra scalar field that represents the volume 
fraction of a specific fluid. This model is also regarded 
as the advection-diffusion- reaction equation dCalder et all 
l2007l : lTownslev et al .1 . 120071). The method is kno wn for its ex¬ 
act conservation of mass ( Scadovelli fe Zaleskil . If 9991 ) . How¬ 
ever, the implicit nature of geometric quantities becomes its 
disadvantage. The geometry of the flame, which is important 
in reconstructing the thermodynamics of partiall y burnt 
grids , is dependent of the reconstruction scheme (iRudmanl 
1 19971 ). In SNIa conte xt, t he algo rithm is pioneered in sim¬ 
ulations presented in iKhokhlovl (119931) . The FLASH code 
has adopted this algorith m as the default flame capturing 
scheme; see for example (ICalder et all [200711 . 


1.4 Sub-grid Turbulence 

The fluid motion can be turbulent because of the high 
Reynolds number. The flame surface is known to be un¬ 
stable subject to hydrodynamic instabilities related to tur¬ 
bulence. The eddies stretch and perturb the flame front, 
which in crease its effect ive burning area and enhance burn¬ 
ing rate <T immesl . ri992l L when compared with laminar flame 
at the same density. Similar to the difficulties in resolving 
the flame interface, the eddies have sizes do wn to the Kol¬ 
mogorov’s scale (~ cm) (IWoosIev et al.l . [20091 1. and therefore 
a direct numerical modeling of these eddies is not yet feasi¬ 
ble. Instead, sub-grid turbulence is employed to mimic the 
small-scale eddies by scaling kinematics from large scales. 

Various implementations of sub-grid turbulence models 
are proposed. They can be classified mainly into the one- 


equation model, two-equation model and Reynolds stress 
model. In all these models, the production/dissipation rates 
are closed by direct numerical simulations, statistical closure 
or dynamical modeling. 

_The one-equation model, first proposed bv lSmagorinskvl 

(119631) . introduces the turbulence kinetic energy q. It first 
ap pear s in the SNIa simulation in iNiemever fc Hillebrandtl 
(I1995ah . which relies o n a statist ical model based on Kol¬ 
mogorov’s scaling llClemend . Il993l ). The model gives break- 
throughs b y pr ovidin g a healthy explosion, as shown in 
iReinecke et alj (l2002al ). However, it is argued that the dis¬ 
sipation m ight n ot be foll owin g Kolmogorov’s one instanta¬ 
neously J Schmidt etldl [200i|). The model is later extended 
in I Schmidt et al. ( 20051 ). which generalizes the eddy produc¬ 
tion or dissipation terms dynamically. 

The two-equa t ion model, first proposed in 
lLaunder fe Spalding! (11974 . introduces two extra quanti¬ 
ties, which model the decay of eddies by an extra dissipation 
term e (q — e model) or a vorticity term uj (q — uj model). 
Both models aim at describing the dissipation of eddies. 
However, the original equations do not guarantee the 
positivity_of q and e (non-realizable). We refer the reader 
for modifications of the equations 


to IShih et all 

that can guarantee the realizability condition. Unlike the 
one-equation model, this model is not applied in common 
SNIa simulations. 

The Reynolds stress model provi des algebraic clo sure to 
the sub-grid velocity correlations (IShih et all Il995l ). How¬ 
ever, not all the variables can be closed by direct numerical 
simulations due to the larger number of coefficients. A sim¬ 
ilar problem can already be found in some extensions of 
the two-equation model. For instance, in the three-equation 
model, an ad hoc timescale is neede d in order to close the 
equation set (lYoshizawa et alll2012l) . 


1.5 Motivation and Structure 

In SNIa literature several codes are commonly used for 
studying the explosion phase. The first one is LEAFS, which 
is known as PROMETHEUS i n earlier publications (see for 
example (IReinecke et all 1 19991) ). This code makes use of the 
piece wise-parabolic method (PPM) (ICollela fe Woodward! 
1 1984) with moving meshes that follow the WD expansion 
(see lRoepke fe Hillebrandtl j2005l ); lRoepkd (120051) for the de¬ 
velopment and early applications). The code models the tur¬ 
bulent nuclear flame using a sub-grid turbulence algorithm 
with the level-set method as flame-tracking scheme. This 
code has been used in large-scale massivel y parallel runs in 
PTD dFink et al.l . l20l4 and DDT models (ISeitenzahl et all 
I 2 OI 2 I). Another hyd rodynamics code is the FLASH code 
( Frvxell et aill2000l) . This code is embedded with adaptive- 
mesh refinement such that the flame structure can be traced 
with fine details. Recently, FLASH has been developed into 
a part of a pipeline in SNIa numerical tool in constructing, 
analyzing an d com paring numerical data with observations 
jLong et all I20l4 . The code uses the advection-diffusion- 
reaction scheme in describing the ash distribution and flame 
geometry but without sub-grid turbule nce. Recen t large- 
scale computation s in the PTD mode ls dLong et all l20l4 
and GCD models (IJorden et al.l . 120121) have used this open- 
source code. We compare our code with LEAFS and FLASH 
in Table [l] 
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There are two motivations for us to build our own hy¬ 
drodynamics code. First, at the time when the code was 
being developed, there was not yet a published SNIa code 
which makes use of a 5th order scheme in spatial discretiza¬ 
tion. Second, we shall apply the hydrodynamics code to var¬ 
ious astrophysical contexts. One of our work in progress 
is to st udy th e effec ts of dark matter on SNIa explosions 
(jLeung et all l2015bi ). By building our own code, it is eas¬ 
ier for us to implement and analyze any additional physics 
component in SNIa simulations. 

In this article, we report numerical tests of our hydro¬ 
dynamics code and comparisons with other benchmark tests 
reported in the literature. In Section [5J we describe the hy¬ 
drodynamics equations and the numerical techniques that 
we implement in our code. In Section[3]we describe the input 
physics for modeling SNIa. We also introduce the radiative 
transfer code by describing its numerical implementation. 
Then in Section [I] we outline the post-processing tools we 
have used in extracting observables from the hydrodynamics 
results, including the tracer particle method and the one¬ 
dimensional moment-integrated radiative transfer code. In 
Section [5] we report our results. We describe the configura¬ 
tions and results of some benchmark one- and two- dimen¬ 
sional tests. We also show two test runs using the PTD and 
DDT models, and we compare our results with some similar 
runs in the literature. In Section [6] we describe the results 
after post-processing, including the detailed nucleosynthe¬ 
sis yields, bolometric light curves and the synthetic spectra. 
In Section [7] we summarize and briefly outline our future 
plans. 


2 METHODS 


2.1 Hydrodynamics 


We have developed a two-dimensional Eulerian hydrody- 
namical code in cylindrical coordinates to model SNIa. The 
equations include 


dpv r 

~dT 

dpVz 

dt 


dp dpv r dpv z 
dt dr dz 


d(pV? +P ) dpVrVz 


dr 


dz 


dpVrVz d{pvl + p) 


dr dz 

dr d(r + p)v r 
dt dr 

(r+p)v r ( 9$ S$\ ^ 

(“■*+»■ aTl+«- 


pVr 
r 

d$ 
°~dL’ 
d$ 
r dz 

d(r + p)v z _ 
dz 


PVr 

r 

pVrVz 


( 1 ) 

(2) 

(3) 

(4) 


In the above equations, (r, z) are the distances from the 
polar axis and the symmetry plane, p, v T , v z , p and r 
are the mass density, velocities in r and z directions, pres¬ 
sure and total energy density of the fluid. The total en¬ 
ergy density includes both the thermal and kinetic parts 
r = pe + |pu 2 , where e is the specific internal energy. <f>, 
the gravitational potential, satisfies the Poisson Equation 
V 2< f> = 47rGp. Q = Qnuc — Q v — Qturb is the heat source 
and loss from nuclear fusion, neutrino emission and sub¬ 
grid turbulence. Note that we deliberately arrange the Eu¬ 
ler equations into this conservative form in order to couple 
them with the Weighted Essential Non-Oscillatory scheme 


(WENO) directly. The terms on the right hand side are re¬ 
garded as source terms. 

We use a high-resolut i on sh ock-capturing scheme, 
WENO dBarth fe PeconinckL 11999i '). for spatial discretiza¬ 
tion. This is a fifth-order scheme, which processes piecewise 
smooth functions with discontinuities in order to simulate 
the flux across grid cells with high precision, while avoiding 
spurious oscillations around the shock. In the code, we use 
a finite-difference formulation of the WENO scheme while 
the hydro quantities in the simulations, such as p, v, e and 
so on, represent the averaged values of the actual profiles 
across the Eulerian grid. For the const ruction of numerica l 
fluxes, we use the WENO-Roe scheme d Jiang fe Shul[l996h . 
First we separate the hydro flux into a positive flux and a 
negative flux by a global Lax-Friedrichs splitting. Then we 
apply the WENO reconstruction on the two fluxes by an up- 
winding scheme. The numerical flux is obtained by summing 
the reconstructed positive and negative fluxes. 

For the discretization in time, we use the five-stage, 
third-order, n on-strong-stability-preservi ng explicit Runge 
Kutta scheme dBarth fc PeconinckL 1199^ 1. Note that in the 
original prescription the three-stage, third-order, strong- 
stability-preserving explicit Runge-Kutta scheme is sug¬ 
gested. However, comparison of these two schemes shows 
that th e latter one has no sta bility advantage and is less ef¬ 
ficient I Wang fc Spiteril . t2007i ~). Therefore, we adopt the new 
time-discretization scheme. 

The gravitational potential is obtained by solving the 
Poission equation V 2< f> = 4nGp. We use the successive over¬ 
relaxation method with the Gauss-Seidel method to solve 
the Poisson equation in cylindrical coordinates. The itera¬ 
tion is stopped by default when the relative changes of all 
grids reach below 1CP 8 . The inner boundaries are assumed 
to be reflecting and the outer boundaries are fixed by us¬ 
ing the Roche approximation, which is given by <f>(r, z) = 
—GM/\Jr 2 + z 2 . We remark that the Roche approximation 
assumes point-like mass everywhere, but in general higher 
multipole terms are needed. In our case, since the outer 
boundaries lie about a few times of the initial stellar ra¬ 
dius away from the center of the star, the higher multi¬ 
poles’ contributions are negligible. We find that the quadru¬ 
ple term gives a less than 1% correction to the monopole 
term throughout the simulations. 


2.2 Level-set method 

We use the level-set method to trace the geometry of the 
flame front. This method is widely applied to other ar¬ 
eas which require a detailed description of surfaces. It has 
b een applied i n SNIa modeling as prescribed in details in 
iReinecke et al.l (119991) . For the sake of completeness, we out¬ 
line the governing equations and its procedure. 

The level-set method introduces an extra scalar field S 
imposed on an Eulerian grid, which satisfies | VS'I = 1. The 
flame front is defined by the zero-contour of the held. We 
also define grids with positive (negative) values of S to be 
full of ash (fuel). Geometrically speaking, the value is the 
minimum separation of a given point to the zero-contour 
surface. 

The flame front propagates by both fluid advection and 
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Table 1. Comparison of the numerical components and structure of our code with other commonly used hydrodynamics codes. 


Scheme 

Our code 

LEAFS (PROMETHEUS) 

FLASH 

Spatial discretization 

WENO 5 th order 

PPM 

PPM/WENO 

Accuracy 

5 th order 

3 rd order 

3 rd _ gth or( J er 

Time discretization 

Runge-Kutta 5-Step 

Directional splitting 

Directional splitting 

Accuracy 

3 rd order 

2 nd order 

2 nd order 

Adaptive Mesh Refinement 

No 

No 

Yes 

Dimensions 

1-2 

1-3 

1-3 

Flame capturing scheme 

Level-set ( default ) 

Level-set 

ADR 


Point-set 



Sub-grid turbulence 

Yes 

Yes 

No 

Number of online isotopes 

19 

5 

3 

Hardware acceleration 

GPU 

Multi-processors 

Multi-processors 


flame burning, which can be written as 

_ = -(« + <% ame )-VS. (5) 

v is the velocity of the fluid, governed by the compress¬ 
ible Euler equations, Ffl ame is the effective flame propagation 
speed obtained from the nuclear reactions described in Sec. 
WI\ which points towards the outward normal direction of 
the flame surface. For the dependence on density or compo¬ 
sition of the flame speed, see in Appendix|A]their derivations 
and analytic fits. 

With the evolution equations, we may simultaneously 
evolve the flame front coupled with the hydrodynamics. Note 
that in our simulations, the scalar field is defined on the 
corners of a grid, such that the position of flame front on 
the grid boundary can be found exactly. 

In updating the scalar field, we make use of operator 
splitting to deal with the contributions from fluid advection 
and flame propagation. The latter part can be done with a 
common high-order advecti on scheme such as the p iecew ise- 
parabolic method (PPM) (jCollela fc Woodward! . 11984 ) or 
the WENO scheme, such that the advection may preserve 
the area (volume) of the flame. Naturally, we use the WENO 
scheme as we have used to solve the hydrodynamics. After 
each update in S, we need to do a reinitialization (or re¬ 
distancing) such that the scalar field S can maintain the 
distance property |VSj = 1. We describe the procedure be¬ 
low. 

First, we enumerate all positions where the flame front 
cut the grid, by finding all neighboring pairs of the scalar 
field where its values are of opposite signs. The position of 
the flame front is obtained by linear interpolation. Then the 
distances of all grid points to the flame front positions are 
compared, and the smallest one is chosen as the new value 
of S. 

Note that this procedure preserves the distance prop¬ 
erty of S, but it perturbs the current flame position as well. 
Therefore, a smooth transition from its original value to the 
corrected value is used. When the grid is far away from the 
flame front, the distance value is used, otherwise the orig¬ 
inal value is kept. Mathematically, for a given scalar field 
S 0 \d(i,j) with a minimum flame distance duj), we have 

^new(ij') — ^(d(i 1 j))*S , old(i,_ 7 ') T ( 1 ^ ) )^(*o) i (6) 


where H(d) is a smooth function satisfying H( 0) = 1 
and H(o o) = 0. In our case, we follow the choice of 
iReinecke et alj (119991 ) that 


H(d) 


1 - tanh W73 

1 - 


(7) 


with do = 3A, A being the grid size in our simulations. This 
formula implies that the values of S on the three closest 
neighboring grids around the flame front are unchanged. 

As we have mentioned, the level-set method is not the 
only flame-capturing method. We have also developed the 
point-set method in order to compare their performances. 
See Appendix [C] for the details and numerical results. 


2.3 GPU accelerations 

In view of a large number of subroutines involved in each 
time step, we share parts of the computation with a graph¬ 
ics processing unit (GPU). By allowing the computations to 
be simultaneously done by more than thousands of threads, 
the total computation time can be drastically reduced, es¬ 
pecially when the job requires repetitive but similar opera¬ 
tions. 

One major difference in CPU coding and GPU coding 
is the lack of static memory in GPU. Such a feature has 
made the use of GPU in typical hydrodynamics simulations 
unfavorable because the static memory provides convenient 
access to data, including the hydrodynamics variables and 
chemical composition, without the need of passing a large 
amount of variables among subroutines. Furthermore, such 
data migration, if not handled carefully, can take up even 
more time than using CPU solely. Nevertheless, the use of 
GPU can be advantageous when applied to subroutines that 
do not have the above undesirable properties. 

In Table [2] we tabulate the subroutines being processed 
by an NVIDIA GT Titan Black GPU and their running 
time, compared with an Intel 17 quad-core CPU. They in¬ 
clude the Poisson solver, reinitialization procedure in the 
level-set method, WENO scheme and spatial discretization 
in other subroutines. 

In general, two classes of time reduction can be found. 
For subroutines with pure matrix operations or repetitive 
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Table 2. A sample run of the code and a comparison of running 
time of subroutines with GPU version in one step. Subroutine 
time lapses are in unit of seconds and the full run is in units of 
hours. CPU: Intel 17, GPU: NVIDIA Titan Black. 


subroutines 

run time 
CPU only 

run time 

CPU + GPU 

WENO scheme (s) 

3.45 

2.06 

Poisson solver (s) 

468.42 

12.07 

Reinitialization (s) 

20.18 

1.20 

Spatial discretization (s) 

1.08 

1.10 

One full run (hour) 

199.66 

29.78 


procedures, such as the Poisson solver or the reinitialization 
scheme, a factor of 20-40 reduction can be obtained. For sub¬ 
routines with heavy data transfer or significant amount of 
logical operations, such as the WENO scheme or the spatial 
discretization (which includes flux limiter), the time reduc¬ 
tion is less significant, about a factor of 2 or less. Overall, 
there is a factor of 6 reduction in run time by using both 
CPU and GPU, compared with runs which make use of CPU 
only. 

3 INPUT PHYSICS 

3.1 Equation of states 

To close the equation set, we use the ope n-sourc e H elmholtz 
equat ion o f state (EOS) repor ted in iTimmes fc Swestvl 
ll 19 thill and ITimmes et al'l (120091 1 . The EOS describes the 
equilibrium thermodynamical properties of a gas, which in¬ 
cludes 1 . electrons in the form of ideal gas with arbitrar¬ 
ily degenerate and relativistic levels, 2 . ions in the form of 
classical ideal gas, 3. photons in Planck distribution, and 4. 
contributions from electron-positron pairs. The subroutine 
takes input of p, temperature T, mean atomic mass A, mean 
atomic number Z and gives other thermodynamical quan¬ 
tities, including p, e, adiabatic index 7 and so on, together 
with their derivatives with respect to each input parameter. 

3.2 Nuclear Reaction Network 

To calculate the nuclear reaction heat production Q nnc , we 
use the 19- isotope nuclea r react ion network subroutine re¬ 
ported in IF. X fc Arnettl dl999l l. The isotopes include 1 H, 
3 He, 4 He, 12 C, 14 N, le O, 20 Ne, 24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 
44 T, 48 Cr, 52 Fe, 54 Fe, 56 Ni, neutron and proton. The fusion 
network includes series of reactions from hydrogen burning 
up to silicon burning. Both (a, 7 ) and (a, p) (p, 7 ) reactions 
are included. Other isotopes, including 2 ' Al, 31 P, 35 C1 and 
so on are included implicitly. Certain isotopes are less impor¬ 
tant here, for example, 4 H, 3 He and 14 N. They are involved 
in hydrogen burning or CNO-cycle, which are not crucial 
in SNIa. On the other hand, isotopes along the a-particle 
chain are important for the exothermic processes in the ex¬ 
plosions, as most energy is released by these processes, and 
the amount of 56 Ni is important for observations. 

We do not couple the nuclear fusion subroutine directly 


in the code except at low density for two reasons. First, 
the deflagration and detonation wave fronts at high density 
have widths much thinner than the numerical grid size, and 
so only a fractional change occurs in the grids which are 
partially burnt in each time step. The energy release and 
isotope variations will be over-estimated if the whole cell is 
considered. Hence, the energy release relies on the fractional 
changes of area (volume) being consumed by the deflagra¬ 
tion and detonation wave front, which depends strongly on 
the local density, but weakly on the local temperature. On 
the other hand, at low density, the deflagration and deto¬ 
nation wave can have a width comparable with simulation 
grid size. Then computing the nuclear reaction yields with 
data from the whole cell is acceptable. Second, the accuracy 
cannot be much improved unless the isotope transport of 
the burning grids are well described. But an exact recon¬ 
struction is time-consu ming and sensit i ve to the numerical 
accuracy. As shown in iReinecke et alj (I1999T ), an exact re¬ 
construction is difficult for grids with large numerical errors. 
As a result, we choose the less accurate scheme by using ta¬ 
bles of reaction products and energy release as functions of 
the input density. We describe the construction method and 
results in Appendix m and Appendix [Bl 

3.3 Sub-grid turbulence 

The use of the Euler equations is a good approximation when 
the fluid viscosity is small. This is true for the case of a WD. 
In a WD the fluid has a typical Reynolds number R e ~ 10 14 . 
This gives a dissipation range 7 related to the integral range 
L by the Reynolds number 

L = V Rl /4 . ( 8 ) 

Therefore, the dissipation range is about lO -3 cm, and the 
scales in which turbulence needs to be directly modeled are 
far smaller than the simulation grid size A « 10® cm. In 
this range p < A < L, called the inertial range, the hydro¬ 
dynamics is dominated by the inertial term in the Navier- 
Stoke Equation and the turbulent dissipation is independent 
of small-scale viscosity. This allows us to treat the fluid as 
an ideal fluid, and take the effects of unresolved turbulence 
statistically. In particular, the turbulent velocity satisfies the 
Kolmogorov’s scaling law, which suggests that the local ve¬ 
locity fluctuations v(l) are related to the local length scale 
l by 

«(Z) =«(£)(£) 1/3 , (9) 

and there is a constant energy transfer rate from large to 
small scales. 

The sub-grid turbulen ce m odel i s bas ed on the above 
properties as described in IClementl (| 1993l l. This model is 
originally used to resolve the problem of frozen flow as com¬ 
monly found in simulations of astrophysical convection. Such 
flow consists of unphysically large scale flow in opposite di¬ 
rections within neighbouring cells. Note that this result is 
also a solution of the Euler equations. It appears because the 
simulation fails to recognize that such large velocity gradient 
can actually result in sub-grid turbulence, which dissipates 
the excess velocity gradient into thermal energy by its kine¬ 
matic viscosity. In practice, the effects of sub-grid turbulence 
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are included as a heat source Q tur b 

• dlJi 

Qturb — —ApqV * V + Yjij — - p^dis ~\~ C^ArchPP'eff ? (10) 

OXj 

where 

E « = '’""(S; + S-- 4V ' rf «) < n > 

is the shear-stress tensor of the fluid flow. CArch is cho¬ 
sen such that the flame propagation speed recovers the 
limit when gravity is important and Rayleigh-Taylor insta¬ 
bilities become the major flame acceleration mechanism. 
We have CArch ~ 1/2 in our case. A is a dimension¬ 
less const ant with A = I (2/3) for 2- ( 3-) dimensional 
fluid f low (Nie mever fe Hillebrandll Il995al : iReinecke et all 
l2002al l . The terms on the right hand side of Eq. m are re¬ 
garded as sources, which include turbulence production by 
compression and shear, turbulence dissipation and Archi- 
median production. g e e is the effective gravitational force 
defined by the product of the local gravitational force and 
the Atwood number. The Atwood number characterizes the 
density contrast between burnt and un burnt matter, w hich 
are functions of density, as presented in lTimmes fc Wooslevl 
(|l992h . q is the specific turbulence kinetic energy q = v 2 / 2 , 
which evolves as a scalar field and couples to the fluid 
through Qturb by 

- + V • (pvq) = Qturb + V • (prWbVg), (12) 

where the last term on the right hand side corresponds to 
turbulence diffusion. 

In this paper, in order to compare consistently with re¬ 
sults from other hy drodynamics codes, we use the algebraic 
closure proposed in lKhokhlov et all dl995l ). The turbulence 
generation/dissipation terms are derived from Kolmogorov’s 
scaling relation, namely 

r-turb « Av(A) = CAq 1/2 , (13) 


and 


Cdii 


A A 


(14) 


with C and D as parameters. In dClementl . Il993l ) C and D 
are modeled in analogy to the ” wall proximity functions” by 
defining the parameter W = t/q. It is found that 


C = 0.1A, 


(15) 


D = 0.5 /F, (16) 

with 

F = min[100, max(0.1, 10 _ 4 W)]. (17) 

The term F ensures that the generation term is large when 
q is small (laminar flow) and the dissipation term is large 
when q is large (turbulent flow). We remark that even we 
have mentioned Kolmogorov scaling at the beginning of this 
section, this sub-grid turbulence model does not assume any 
particular scaling relation, since the forms of Eqs. CUB and 
Cl) are obtained by only dimensional analysis. Individual 
scaling relations in the fluid flow can be realize d by c arefully 
tuning Eqs. (fl5l) - (fTTll . But as indicated by [clement] d 19931 ). 
this procedure can be very difficult unless one solves the 
spectral equations simultaneously, which is unfeasible in the 


current stage. Nevertheless, the current form is shown in the 
literature that it can describe astrophysical flow with a high 

Reynolds number well. _ _ 

Following the suggestion of lSchmidt et all (2006), the 
turbulent flame propagation speed is connected to the tur¬ 
bulent velocity fluctuation by 


Aflame — ^lam 



(18) 


in the asymptotic regime of turbulent burning dPocheaul . 
119941) . with Ct = 4/3. It is easy to observe that at q = 0, 
Afl ame = flam automatically and at g » v ? arr , , Uflame ->• 
\f8q/3. The laminar flame speed is obtained by solving 
for the structure of a steady state deflagration wave. The 
method and numerical fit are presented in Appendix m 


3.4 Neutrino Emission 

We use a open-source thermal neutrino emission subrou¬ 
ting which calculates the neutrino luminosity Q„ at a given 
temperature, density, mean atomic number and mean pro- 
ton number, using the analytic fit presented in lltoh et aid 
d 19961) . Major thermal neutrino production channels are in¬ 
cluded, such as the pair-, photo-, plasma-, bremmstrahlung 
and recombination neutrino processes. 

To construct the neutrino spectra, we consider the neu¬ 
trino emissivity due to different neutrino generation mech¬ 
anisms A neu = A pa , ir + Apiasma + ■•• ■ In our calculation, 
we focus on two major mechanisms: pair-annihilation pro¬ 
cess and plasma neutrinos. The pair-annihilation neutrino 
emissivity a t relativistic and non-d egenerate regimes has an 
analytic fit (Misiasz ek et. aid . 120061 ) given by 

Apair. = F</>(E V ), (19) 

with 

F = G Z m f (Ml 0 + Ml 0 ) (20) 

lo7r & 

being the total number emissivity and 

<K E *) = ^ (H) exp (~aE„/kT) (21) 

being the shape function. E v is the neutrino energy. m e is 
the electron mass in MeV. Ai, a and a are obtained from 
parametric fits of the exact relations. Gf is the Fermi weak- 
coupling constant, while Ml 0 ( M +°) is the zeroth moment 
of electron (positron) energy defined by 


M nm 

= (7C 2 v - 2C\)Gf /2 _ 1/2 Gl /2 _ 1/2 + 



9C'yG*/ 2 G ^/ 2 + (Cy + C\) x 


(4^71/2+1/2 Cm/2+1/2 ^n/2-1/2^771/2+1/2 



^n/2+1/2 ^ 771 / 2 - 1 / 2 )’ 

( 22 ) 

with 

G n (a,/3) a 3+2n J a 1 + exp (x±/3) dX 

(23) 


1 available in http://cococubed.asu.edu/code_pageseos.shtml 
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being the Fermi integral and a, /3 and x are the dimension¬ 
less electron mass, chemical potential and energy. Cv and 
Ca are the vector and axial coupling constants. 

There is a l so an analytic fit for plasma neutrinos given 
in lOdrzvwolekl (|2007l ) 

N plasma (^l') = A 2 kTmtexp(-E„ / kT), (24) 


where 


s~l2 

8n 4 ah 4 c 9 


(25) 


and mt is the transverse electron mass related to the Fermi- 
Dirac distributions of electron fi and positron / 2 by 


2 4 a 

m t = — 

7T 


F 


E 


(/i + h)dp. 


(26) 


In our simulations, the transverse electron mass and 
the pair-annihilation total neutrino emissivity are tabulated 
as functions of density and temperature. The emissivities 
of these two channels can then be readily computed and 
summed to obtain an energy distribution of neutrinos. 


4 DATA POST-PROCESSING 

This section introduces the data processing after the hy¬ 
drodynamics runs. They include the tracer particle algo¬ 
rithm for reconstructing detailed nucleosynthesis yields, and 
a moment-integrated radiative transfer code for predicting 
the corresponding bolometric light curves. 


4.1 Particle Tracer 

The use of the 19-isotope network reaction products in tab¬ 
ular form only provides an approximation because of the 
absence of electron capture and other off-a-chain isotopes. 
In the literature, table forms of the energy production and 
chemical composition are often used due to the significant 
difference between reaction and hydrodynamics time scales 
and length scales. However, such results are inadequate if 
we need to compare the nucleosynthesis with real observa¬ 
tional data, which can be recorded in fine details. In order to 
obtain a precise distribution o f elements, we us e the t racer 
particle method prese nted in ISeitenzahl et alJ (I2004T ) and 
ISeitenzahl et al.1 (|20ld 'l. This method introduces a number 
of pseudo particles of either equal mass or equal volume. 
They follow but not affect the fluid flow. The density and 
temperature of these particles are recorded as functions of 
time. In simulations, the particle density, temperature and 
velocity are obtained by bilinear interpolation from the fluid 
properties. 

After the hydrodynamics simulations, the density- 
temperature trajectory of each particle is recalled to trace 
back its chemical evolution. In our simulation, the recon¬ 
struction is done based on the same 19-isotope network, 
which can be compared with the results obtained from hy¬ 
drodynamics runs. 


4.2 Radiative Transfer for Bolometric Light Curve 

After the simulations, we map the hydro results from the 
two-dimensional Eulerian form to a one-dimensional La- 


grangian form with spherical symmetry. During this trans¬ 
formation, the total mass, energy and momentum are con¬ 
served in order not to provide spur ious perturbat ion to 
the pr ofiles. We follow the method in IZhang fe Sutherland! 
dl994h to construct the light curves by solving the one¬ 
dimensional time-dependent moment-integrated radiative 
transfer equations, which include two equations from La- 
grangian hydrodyanmics, 


Dt \p) or c r z 

+ (p + J'ar)-^- 0^ = CK E E r - 4nnpB(T) + Q decay, (28) 
and two equations from radiative transfer, 



+ 


D 


— ( 3 / — 1 )- 


J Dt 

A-kkpB — ckeEt — 


E r ~ 


pr 

d(4ivr 2 Fr 

dm 


(29) 

(30) 


1 DFr 1 d{fqE r ) 
c 2 Dt + q sph dr 



( 31 ) 


The hydrodynamics part is almost identical to those in 
previous parts. However, the physical quantities are defined 
in a staggered grid as typical in Lagrangian hydrodynam¬ 
ics. Density p, specific internal energy e, fluid pressure P, 
enclosed mass m(r), specific energy production rate due to 
the decay of radioactive isotopes Qdecay, radiation-energy- 
mean opacity k e , Planck-mean opacity Kp, flux-mean opac¬ 
ity xf, blackbody radiation rate B, Eddington factor / and 
sphericity q sp h are defined on grid centers (r = r„_|_i/ 2 , n = 
1, 2, 3...). Velocity v and radiation flux F r are defined on 
grid boundaries (r = r n , n = 1, 2, 3 ...). 

iz a r is the artificial viscosity defined on the grid cen¬ 
ters, which operates whenever fluid compression occurs. At 
grid centers k + 1/2, zz ar is proportional to the difference in 
the magnitudes of the velocities of nearest grid boundaries, 
namely 


v a,T k+l/2 


= 2 


/ n n\2 n 

-i( w fc + l — v k ) Pk+ l/ 2 i 


(32) 


whenever D(l/p)/Dt < 0 and dv/dr < 0 and i/ num = 4 is the 
numerical viscosity coefficient, which controls the smearing 
of shocks. The blackbody radiation rate is given by Planck’s 
formula B(T) = (ac/'in)T 4 . ke and Kp are the radiation- 
energy-mean opacity and Planck-mean opacity, given by 


j K(v)Evdv/ 1 E u dv, 

(33) 

Kp = J K(v)B v (T)dv, 

(34) 


with B„(T) = 2 his 3 /(e x — 1) being the Planck distribution 
function at some given temperature T, frequency v and x = 
hu/kT. \F is the radiation-flux-mean opacity, given by 


Xf 


/ 


X(i')F v dv/F r . 


(35) 


Since we consider only bolometric luminosity, which means 
all physic s quantities are freque ncy-integrated, we follow the 
choice of iHoeflich et all (11993h that ke = Kp, and xf is 
replaced by the Rosseland mean opacity Kp, defined by 


kr 


f^du ■ 


(36) 
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We describe how we model the continuum opacity source 
and include the atomic transition lines as opacity sources in 
Appendix El 

The Eddington factor is an algebraic closure of the ra¬ 
diative transfer equations, given by 

fu = T»dpJ / (/ ludn J , (37) 

with fi = cos# being the direction cosine of radiations with 
respect to the radial direction. I v is the monochromatic ra¬ 
diation intensity at frequency v. The sphericity q sp h(i') mea¬ 
sures the level of isotropy of light rays, defined by 

In (qA = f [(3/„ - l)/(r'f u )dr'] , (38) 

J r c 

with r c being the radius, within which the matter is optically 
thick. 

In our simulations, Q decay is the energy release due to 
decays of 56 Ni and 56 Co. In the radiative transfer, we con¬ 
sider only the time evolution of three isotopes: 56 Ni, 56 Co 
and 56 Fe. The energy release is related to the decay rates 
given by 


Q decay — Dry 


DX m 

1 Dt 


+ £ Co 


DXco 

Dt 


, e+ DX Co 

+ £ Co , (39) 


where and £q 0 are the specific energy release due to the 
decays of the isotopes by emitting gamma rays, and is 
that by emitting positrons. The amount of isotopes can be 
computed analytically, 


Am = X Ni(ini) exp(—t/r N i), (40) 
Aco = A N i( ini )——-[exp(— t/rco) - exp(-t/rNi)], (41) 

T’Co TNi 

with TNi = 7.605 x 10 s s and rco = 9.822 x 10 6 s the de¬ 
cay half-lives of the two radioactive isotopes. Am and Aco 
are the mass fractions of nickel-56 and cobalt-56. The energy 
released by the positron decay channel is assumed to be com¬ 
pletely absorbed by local matter because a positron has a 
much higher scattering cross section than a photon, even at 
low density. Only a fraction of photons is absorbed, which 
is controlled by the deposition function D ' 7 . To determine 
the d epositi on fu nction , we follow the prescription described 
in [Swartz et alj (jl995f ) by solving analytically the grey ra¬ 
diative transfer of the gamma ray radiations in the two- 
stream approximations. Specifically, we obtain the energy- 
integrated intensity for both incoming (/“) and outgoing 
directions ( I + ) by solving 
± 

±-^ = v -K,pI ± , (42) 


subject to boundary conditions I~ = 0 at the surface of the 
ejecta and I~ = I + at the core, r] is the frequency integrated 
gamma ray emissivity. After solving for 7* at each grid shell, 
the deposition function is computed by 


Dj 


AnK-yJ 

Q decay 


(43) 


with J being the moment-integrated intensity 



(44) 


which is integrated over all solid angles dtt. 


In Eq. (1371) . we need to solve /„ formally. We decom¬ 
pose it into the symmetric part j v = [7„(r, p) + 7„(r, — p)\/2 
and anti-symmetric part h v = [Iv(r,p) — Iv{r, —p)]/2. They 
satisfy 


1 Dj v dh v 
c Dt ds 


KvpB v (T) — 


K v p+ (1 — 3 p~) - 

rc 


(1 + p 2 ) Dhip 


Dt 


1 Dhv dj v 
c Dt ds 


-Xvphv 



(45) 

(46) 


In the above equations, s 2 = r 2 — p 2 where p is the impact 
factor of the ray. By solving Eqs. (145II and (1461) with bound¬ 
ary conditions j v (r = R) = h v (r = R) and hv(r = 0) = 0 
for all s, we may solve for /^(r, u), which then provides the 
Eddington factor f v and sphericity q v . 

To solve the radiative transfer equations, we input 
the density, velocity, temperature and composition profiles. 
Since the explosion does not yet achieve the homologous 
expansion, we map the initial velocity profile to a homolo¬ 
gously expanding one with the same kinetic energy by using 
the results from the end of the simulations of the explosion 
phase. Certainly, a consistent way to obtain a homologous 
profile is to let the system evolve. This requires either a suf¬ 
ficiently large simulation box or box size that varies with 
time. However, the first way requires an impractically large 
amount of computational resource, while most of it is not 
used except at later time when the star starts to expand. On 
the other hand, with the second way we can keep track of the 
fluid motion while using a manageable computer resource, 
but it requires special numerical treatment of physics com¬ 
ponents that are sensitive to the resolution, including the 
sub-grid turbulence and the level-set method. Therefore, at 
this stage, we artificially replace the velocity profile by a 
homologous one which conserves the total energy of the sys¬ 
tem. 


5 HYDRODYNAMICS RESULTS 

This section presents the hydrodynamics results of var¬ 
ious code tests. They include one-dimensional and two- 
dimensional code tests, such as the standard shocktube tests, 
Gresho vortex and tests of hydrodyanmics instabilities. They 
aim at testing the validity of the code and its capability in 
shock-capturing, level of accuracy and numerical diffusion. 
Then we present hydrodynamics from explosion models in¬ 
cluding the PTD and DDT models. Though it is known that 
the PTD model cannot provide a healthy explosion that 
matches typical SNIa ejecta, its slowness in flame propa¬ 
gation allows various hydrodynamics instabilities to form, 
which can be used to check the physical components collec¬ 
tively, such as the level-set method, sub-grid turbulence and 
products of the deflagration wave. The DDT model is one 
of the possible SNIa mechanisms. Also, our results can be 
compared directly for models with similar configurations in 
the literature. 


5.1 One-dimensional Code Test 

In this two-dimensional code, we study the one-dimensional 
limit by reducing the number of grids in either one dimen- 
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Table 3. Input parameters for one-dimensional code tests. 


Test 


XQ 

Pl 

VL 

Pl 

PR 

vr 

PR 

1 

0.2 

0.3 

1.0 

0.75 

1.0 

0.125 

0.0 

0.1 

2 

0.5 

0.5 

1.0 

-2.0 

0.4 

1.0 

2.0 

0.4 

3 

0.035 

0.4 

5.99924 

19.5975 

460.894 

5.99242 

-6.19633 

46.0950 

4 

0.012 

0.8 

1.0 

-19.59745 

1000.0 

1.0 

-19.59745 

0.01 




0.5 1 

position 


Figure 1. Density, pressure, velocity and internal energy in Test 
1 at tg = 0.2 and xg = 0.3. 



Figure 2. Density, pressure, velocity and internal energy in Test 
2 at to = 0.5 and xg = 0.5. 


sion to a few grids and assuming Cartesian coordinates, i.e. 
the volume of each grid is independent of the radial dis¬ 
tance from the axis of symmetry. We follow the tests pro¬ 
vided in ( iTord . l2007h which are all shocktube tests to under¬ 
stand whether the code can correctly evolve various types of 
shocks. They all have exact solutions computed by Riemann 
solvers. In this part, we use a polytropic EOS with the ratio 
of specific heat 7 = 1.4. The EOS is written as p = pe( 7 — 1). 
The simulation takes place in a spatial domain of interval 
[0,1] with 2000 computing cells. The Courant number is 
fixed to be 0.5. The initial data consists of two parts, the left- 
hand state ( Pl,vl,Pl ) and the right-hand state (pr,vr,pr). 
The transition takes places at xg and the simulation is run 
for a duration tg. We tabulate the input configurations in 

Tabled 

In Fig. □ we plot the density, pressure, velocity and in¬ 
ternal energy of the fluid for Test 1 at t = 0.2. This test is 
a modified version of Sod’s test which assesses the entropy 
satisfaction property of the numerical scheme. The final so¬ 
lution includes a right shock wave, a right-traveling wave 
and a left sonic rarefaction wave. As seen from the figure, 
the code can very well preserve the sharp edge of the shock 
and contact discontinuity. Also, the entropy glitch, which 
usually appears in low-order schemes in a bump inside the 
rarefaction, does not appear. 

In Fig.[2j we plot similarly the density, pressure, velocity 
and internal energy of the fluid for Test 2 at to = 0.15. 
This test contains a solution of two rarefaction waves and a 
stationary contact wave. The solution contains also a low- 
density region, and it tests how well the code can handle 
low-density fluid flow. The numerical solution follows the 
analytic solution very well, except for an undershoot in the 



position 


Figure 3. Density, pressure, velocity and internal energy in Test 
3 at to = 0.035 and xg = 0.4. 


velocity near x = 0.5 and an overestimation by a factor of 
2 in the internal energy. Despite these, the performance of 
the WENO scheme is in general much better than other 
low-order schemes, which result in spurious oscillations in 
the velocity held and a larger overestimation in the internal 
energy (See ITord d2007h for a detailed comparison among 
different numerical schemes). 

In Fig.[3]we plot the test results for Test 3 at to = 0.035. 
This test assesses the robustness of a code in processing 
strong discontinuities. There are two right-traveling shock 
waves, a right-traveling contact discontinuity, and a left- 
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Figure 4. Density, pressure, velocity and internal energy in Test 
4 at to = 0.012 and xg = 0.8. 

traveling shock. It can be seen from the figure that all dis¬ 
continuities can be resolved without smearing. 

In Fig. [4] we plot the density, pressure, velocity and in¬ 
ternal energy in Test 4 at to = 0.012. This test challenges 
the code in evolving a slowly-moving or stationary shock 
discontinuity. The final solution contains a left-traveling rar¬ 
efaction wave, right-traveling shock wave and a stationary 
contact discontinuity. The numerical solution largely follows 
the analytic solution except for a small overshoot in pres¬ 
sure, velocity and internal energy at x = 0.4 at the edge of 
rarefaction wave. Also, there are spurious oscillations with 
a small amplitude in density at x = 0 . 8 . 

Combining these four tests, we can see that our code 
performs well in treating numerical discontinuity and it can 
mimic closely the analytic solution. Typical numerical inac¬ 
curacies, including the entropy glitches, spurious oscillations 
and smearing are highly suppressed. 

5.2 Two-dimensional Code Test 

The one-dimensional tests aim at validating the shock¬ 
capturing properties of the WENO scheme under different 
situations. In this part, we extend the code test to two di¬ 
mensions. Besides shock-capturing, we also study the nu¬ 
merical diffuseness and the advection of smooth flow of 
the code. In all tests the polytropic EOS identical to the 
one-dimensional code tests is used with specific heat ratio 
7 = 1.4. The Courant numbe r is taken to b e 0.5. 

The first test Gresho ( Ores hot, liflflrif is to evaluate 
the code performance in handling time-independent solu¬ 
tion, which tests the code in handling advection of smooth 
functions. The simulation takes place in a square box of 
—0.5 < x < 0.5 and —0.5 < y < 0.5 with A = 0.02, 
0.01, 0.005 and 0.0025 units. The boundary is free every¬ 
where. The simulation time is 3 (code unit). The initial 
profile is a steady solution given by p = 1 everywhere, 
v(r) = 5r, p(r) = 5 + 25/2r 2 at 0 < r < 0.2, v(r) = 2 — 5r, 
p(r) = 9 + 41n(r/0.2) + 25/2r 2 — 20r at 0.2 < r < 0.4 and 
v(r) = 0 and p(r) = 3 + 41n2 at 0.4 < r < 0.5. Here r is 
the distance from the origin. Obviously, any departure of 
the final solution from the initial one must be due to nu¬ 
merical errors from the advection and numerical diffusion. 
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Figure 5. Density contour plot of Gresho test at t = 3. 



Figure 6. Density contour plot of Blast test at t = 1. 

In Fig. [5] we plot the linear density contour plot at t = 3 
for dt = 0.005. It can be seen that the circular structure is 
maintained. The outer part remains in uniform pressure. In 
Table [4] we list the numerical errors in total mass, energy 
and the L 1 -norm of the fluid pressure. At high resolutions, 
the results show a third order convergence which is consis¬ 
tent with the order of the time-discretization scheme of the 

code. _ _ 

The second test is the Blast test riTord. 120071') which 
takes place in a rectangular box of dimensions —0.5 < x < 
0.5 and —0.75 < y < 0.75 with reflecting boundaries ev¬ 
erywhere. The resolution is 512 x 768. The initial condi¬ 
tion is given by p = 10 at r < 0.1 and p = 1 at r > 0.1. 
p = 1 everywhere. This test studies again the capability of 
the code in handling the low density region, the sharpness 
of the hydrodynamical instabilities and the shock-shock, or 
shock-contact discontinuity interactions. In Fig. QT] we plot 
the linear density contour at t = 1.0. It can be seen that 
the symmetry along the x-axis and along the y-axis are pre- 
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Table 4. Fractional changes in mass A M/M and energy A E/E in Gresho test under different resolutions. (p — p rB f) is the L 1 norm 

of the error in the fluid pressure, with reference to the initial state. 


A 

N x 

Ny 

A M/M 

A E/E 

T^p-Pref) 

2 x 10“ 2 

50 

50 

2.96 x 10“ 3 

4.08 x 10“ 4 

4.50 x 10 -3 

1 x 10 -2 

100 

100 

1.09 x 10 -3 

1.52 x 10 -4 

1.63 x 10 -3 

5 x icr 3 

200 

200 

1.49 x 10“ 5 

2.10 x 10“ 5 

2.89 x 10“ 4 

2.5 x icr 3 

400 

400 

1.94 x 10“ 6 

2.73 x 10“ 6 

3.26 x 10 -5 
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Figure 7. Pressure contour plot of Explosicm test at / = 2.667. 



Figure 8. Density contour plot of KH test at t — 1. 

served. Also, the Richter-Meshkov instability at the bound¬ 
ary of contact-discontinuity in the form of dense filaments 
or fingers are sharply captured. 

The third test is the Explosion test dTorol . 120071) which 
is very similar to SNIa explosion but without nucleosyn¬ 
thesis. The simulation box is 512 x 512 with dimensions 
0 < x < 1.5 and 0 < y < 1.5. The inner boundaries are 
reflecting and the outer boundaries are free. The initial con¬ 


dition is given as an overpressure region at the core, with 
p(r ) = 1.0 and p(r) = 1.0 at r < 0.4 and p(r) = 0.1 and 
p(r) = 0.125 at r > 0.4. In Fig. [3 we plot the linear density 
contours at t = 2.667. At this moment, the weak reflected 
shock has arrived inside the shock sphere, and the unstable 
contact discontinuity of the surface in the form of plumes 
can be observed. Again, a sharp symmetry with respect to 
a reflection about x = y can be observed. 

The last test is the Kelvin-Helmholtz (KH) test 
ll McNally et alJ . 1201 il l which studies the KH instability. The 
simulation box is 256 x 256 with dimensions 0 < x < 1 
and 0 < y < 1 with periodic boundaries everywhere. The 
simulation is done up to t = 2. The initial condition is 
that two fluids of different densities flow in opposite di¬ 
rections with a perturbation: p — 1 everywhere, p = 2 at 
0.25 < y < 0.75 and p = 1 otherwise. An initial velocity is 
given as v y = A( 1 + cos(27rx))(l + cos(27 ry)). In Fig. [ 8 ] we 
plot the linear density contours at t = 1. The KH instability 
in the form of spiral shape outside the denser region can be 
observed. 

Combining these four tests, we have seen that the code 
performs satisfactorily with the desired accuracy and shock¬ 
capturing ability in multi-dimensional runs. Also, the code 
has shown its ability in handling hydrodynamical instabili¬ 
ties, such as the Kelvin-Helmholtz instability. 

5.3 Code Test for PTD models 

In this part, we return our focus to code test of SNIa ex¬ 
plosions using t he PTD models. The sim ulation uses the 
Helmholtz EOS llTimmes &; Swestvl . 1 1 9991 ) and all the sim¬ 
ulations are done in cylindrical coordinates with 500 grids 
in each direction. Each grid has a size A = dr = dz ~ 11 
km. A hydrostatic equilibrium model with a given density is 
constructed with an isothermal profile of T = 10 8 K and a 
composition of X( 12 C) = A'( 16 0) = 0.5. The code enforces 
a minimum density of 10 5 g end 3 , which is also the atmo¬ 
spheric density. The time step is limited by the Courant- 
Friedrich-Lewy condition: 

Af max — c c fiA/(c s + \v r \ + \v z \), (47) 

with c c fi = 0.5 and c s the local sound speed. 

An initial three-finger flame front is imp osed, which is 
compa rabl e with the c3 f l ame r eported in iNiemever et al.l 
d 19961) and lReinecke et all (1 19991 ). The reason we choose this 
flame front is the same as that in the literature: we want to 
bypass the slow laminar flame stage and consider the stage 
where Rayleigh-Taylor instabilities are important. The mat¬ 
ter enclo sed by the flame front is f irst b urnt. We follow the 
choice in INiemever fe Hillebrandtl dl995bl ) to set the initial 
specific turbulence kinetic energy q = 10 10 cm 2 s -2 . We also 
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Table 5. Simulation setup for the test of PTD models. Lengths are in units of km and densities are in units of 10 9 g cm 3 . Energy is 
in units of 10 59 erg and masses are in units of solar mass. 


Model 

dx 

Pc 

mass radius 

ExiUC 

-Fkin 

M Ni 

2D-la-PTD 

11.0 

3.0 

1.377 1475.0 

7.141 

2.109 

0.320 


set this value to be the minimum q so as to avoid unphysical 
result when calculating the production and dissipation of q. 

We tabulate the simulation setups in Table. E We plot 
in the upper panel of Fig.|5]the total energy against time for 
Models 2D-la-PTD. At early time the energy release relies 
on slow laminar flame, thus the energy growth is slow. At 
about t fa 0.5 s, the sub-grid turbulence has developed and it 
boosts the propagation of flame and enhances the consump¬ 
tion of fuel significantly. At about 0.7 s, the nuclear flame 
has released enough energy to balance the gravitational en¬ 
ergy. At t fa 0.8 s, the expansion of the star leads to density 
decrease at the flame front, which lowers the energy out¬ 
put of the flame as well. So, the total energy levels off and 
reaches a constant at about 2 x 10 5 ° erg. Notice that our 
initial model and related physi cs are chosen to mimic those 
of the Model c3 — 2d — 256 in iReinecke et ahl (l2002al f . Our 
results are comparable with theirs, which also give the final 
energy at 2 x 10 5 ° erg. In the lower panel, we plot the total 
turbulence kinetic energy q against time. Similar to the neu¬ 
trino luminosity, the sub-grid turbulence energy reaches its 
maximum at around 0.6 s. We plot in the upper panel of Fig. 
EH the neutrino luminosity against time. The model shows 
a single peak at about 0.5 s. This implies that the burning 
is the most vigorous at that moment. In the lower panel 
we plot the neutrino energy spectra with energy 1 MeV - 
5 MeV as functions of time. In early time, the neutrino en¬ 
ergy peaks at 2 MeV. At later time, when the WD starts 
to expand, the 1 MeV neutrino flux becomes the dominant 
one. 

We also plot in Fig. |TT] the time evolution of chemi¬ 
cal isotopes. At early time, iron-peaked elements are pro¬ 
duced, and a considerable amount of 4 He is produced. At 
later time, when the fluid element expands and its temper¬ 
ature drops, 4 He recombines to 56 Ni again, and incomplete 
burning causes the production of a trace amount of IME. 
The results can also be compared with Model c3 — 2d — 256 
l|Reinecke et al.l . l2002al f which observes a total of 0.109 Mq 
IME and 0.40 Mq iron-peaked elements. We obtain sim¬ 
ilar amount of iron-peaked elements but the IME abun¬ 
dance is lower than theirs. This may be due to the dif¬ 
ference in approximating the ash composition and the en¬ 
ergy release at lower density, which is the main site for the 
production. We plot in Fig. [12] the flame shape at t = 1 
s. The initial flame is the c3 f ront shape as described in 
iNiemever fc Hillebrandtl (|l995bl ). which is a ’’three fingers” 
shape. The flame shows signatures of hydrodynamical insta¬ 
bilities: the Rayleigh Taylor instability in the form of mush¬ 
room shape at the top of the ’’fingers” and KH instability 
in the form of curly shape along the ’’fingers”. Injection of 
fuel into the flame can also be seen between ’’fingers”. 



Figure 9. Upper panel: Total energy against time for Model 2D- 
la-PTD. Lower panel: Total turbulence kinetic energy against 
time for the same model. 



Figure 10. Upper panel: Total neutrino luminosity against time 
for Model 2D-1D-PTD. Lower panel: Neutrino energy spectra 
against time for the same model. 

5.4 Code test for DDT models 
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Figure 11. Isotope mass fraction against time for Model 2D-la- 
PTD. 



Figure 13. Upper panel: Total energy against time for Model 
2D-la-DDT. Lower panel: Total turbulence kinetic energy against 
time for the same model. 
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Figure 12. The deflagration and detonation fronts at t = 1 sec¬ 
ond for Model 2D-la-PTD. 

In this part we present results for the DDT model. The 
DDT model is actually a natural continuation of the PTD 
model for two reasons. First, before the detonation starts, 
the evolution is exactly the same as the PTD model. Sec¬ 
ond, the criterion for the first detonation spot relies on the 
local turbulence that the flame width S equals the Gibson 
length, also known as the turbulence length scale iGib, be¬ 
low which the laminar flame is the dominant propagation 
channel. Mathematically, we have 

S = Ig ib = A(i)i am /2g) 3 / 2 . (48) 

The flame width is obtained from the results by solving 
the deflagration wave, as described in Appendix |XJ which 
is a function of density. We note that this relation, which 
is equivalent to the Karlovitz number Ka = 1, only states 
that turbu lence starts to destroy the flame front. But as 
argued in lNiemever fe Wooslevl d 19971 ). this mechanism can 
allow the heat from the ash to be transported to the fuel 
much more efficiently than only by conduction. This creates 
a much wider region with sufficient temperature for carry¬ 



Figure 14. Upper panel: Total neutrino luminosity against time 
for Model 2D-la-DDT. Lower panel: Neutrino energy spectra 
against time for the same model. 

ing out carbon burning simultaneously, which is believed 
to be one of the keys for triggering the detonation. In the 
simulation, whenever there is a grid on the flame front sat¬ 
isfying this criterion, a spherical detonation spot of radius 
A is artificially placed and evolved similarly as the deflagra¬ 
tion front. However, because the deflagration front has burnt 
partially the stellar material in the core, whenever the det¬ 
onation front reaches the flame surface, we assume that the 
detonation front stops. 

We carry out tests for models using DDT and we tab¬ 
ulate the input parameters and final explosion energetics 
in Tabic [6] Th e explos ion energe t ics c an be compared with 
Model Z2>a in lOolombek fc Niemeveil d2005l ) . But the major 
difference between Model Z3a and ours is that the former 
one allows the detonation front to pass through burnt re¬ 
gions, which is not allowed in ours in view of some direct 
numerical simulations of th e instability of the shock f ront in 
the presence of ash product dMaier fe Niemeveid . l200(I . They 
reported E nuc = 14.4 x 10 50 erg with Skin = 8.70 x 10 50 erg. 
The DDT transition time is 0.85 s and the final amount of 
56 Ni in their model is 0.65 M@. Our model shows a higher 
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Table 6. Simulation setup for the test of DDT models. Lengths are in units of km and densities are in units of 10 9 g cm 3 . Energy is 
in units of 10 so erg and masses are in units of solar mass, ttran is the first DDT time in seconds. 


Model 

dx 

Pc 

mass radius 

Enuc 

Ekin 


^tran 

2D-la-DDT 

11.0 

3.0 

1.377 1475.0 

15.651 

10.618 

0.733 

0.70 



Figure 15. Isotope mass fraction against time for model 2D-la- 
DDT. 


2500 



Figure 16. The deflagration and detonation fronts of Model 2D- 
la-DDT at t = 1.00. The temperature distribution is represented 
in different colours. 


energy release, a higher 56 Ni production and an earlier trig¬ 
ger in DDT. 

In the upper panel of Fig. 1131 we plot the total energy 
against time for Model 2D-la-DDT. Prior to the detonation 
front being triggered, the evolution is identical to that of the 
PTD model. By comparing with Model 2D-la-PTD, we see 
that our model is consistent with models in the literature 
that most energy is released by the detonation front instead 
of the deflagration front, which is also needed to match the 
ejecta velocity with the observed SNIa. In the lower panel, 
we plot the total turbulence kinetic energy. The turbulence 
constantly develops as the fluid motion becomes chaotic, and 


then it drops slightly when the WD expands. Turbulence 
grows again when the detonation is launched, which reaches 
its maximum as the explosion has burnt all the fuel in the 
WD. But it then drops sharply due to the rapid expansion 
of the star. 

In the upper panel of Fig. [14] we plot the neutrino lu¬ 
minosity of Model 2D-la-DDT. The neutrino luminosity is 
decreasing when the detonation is just started. But it goes 
up again due to a large amount of matter being thermal- 
ized. In the lower panel we plot the neutrino energy spectra 
with energy 1 MeV - 5 MeV as functions of time. Again the 
number flux of each energy band shows a rise after DDT has 
started and after t = 0.7 s, the 1 MeV neutrino number flux 
exceeds the other four fluxes of higher energy bands. We 
also plot the isotope abundances against time in Fig. 1151 
The early time evolution is comparable to the PTD model. 
After the detonation has started, much more iron-peaked 
elements and IME are produced, leaving trace amounts of 
12 C and 16 0. At last, we plot deflagration and detonation 
flame front at t = 1 s in Fig. 1161 The flame front is shaped 
by the fluid flow, that ’’mushroom caps” and injection of 
fuel can be observed from the flame front. The detonation 
front first starts at the outermost part of the flame front, 
which develops at first spherically, but then loses its sym¬ 
metry when the front encounters another detonation front 
and the deflagration front. 


6 OBSERVATIONAL PREDICTION 

We have demonstrated in previous sections that our hydro¬ 
dynamics code can provide explosion energetics compara¬ 
ble with models in the literature. To further study the hy¬ 
drodynamics result, we perform a series of post-processing 
procedures to extract SNIa observables, including the nu- 
cleosythesis yields by the tracer particle method and the 
bolometric light curve using a moment-integrated radiative 
transfer code. 

6.1 Nucleosynthesis Products 

We use the tracer particle method to keep track of the den¬ 
sity and temperature of each fluid parcel, which represents 
a fixed amount of mass elements in the fluid. We repeat 
the 2D-la-DDT model but with different numbers of tracer 
particles. We study its numerical convergence by adjusting 
the initial radial and angular separation of the fluid parcels. 
In the first test, we use a simple 19-isotope network with a 
fixed output time of hydrodynamics properties at 0.025 s. 
We tabulate the nucleosynthesis yields in Table [T] We also 
include the chemical composition from the hydrodynamics 
run for comparison. In general, the post-processed yields are 
different from that of the hydrodynamics run in two ways; 
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Figure 17. Bolometric light curves of Model 2D-la-DDT. 

the final distributions of 12 C and 16 O are not the same in the 
post-processed results. This is related to the slow-burning in 
the low density regime, which is not assumed in the deflagra¬ 
tion wave product. This also affects the IME distribution so 
that much higher mass fractions of 28 Si and 32 S are observed. 
Second, the amount of iron-peaked elements are drastically 
different in that the 54 Fe are generally lower, showing that 
lower amount of matter has reached nuclear statistical equi¬ 
librium for a complete burning comparing to the expectation 
from the deflagration wave solution. We also compare how 
the number of tracers along each direction affects the final 
results. We find that for N = 80 2 the changes in the final 
yields differ by less than 1 %, while more abundant isotopes 
such as 54 Fe and 56 Ni generally converge faster than less 
abundant elements such as 40 Ca and 44 Ti. 


6.2 Bolometric Light Curve 

Using the hydrodynamics results from models as we have ob¬ 
tained in Sections EH and ED we can predict the expected 
observational signals of these explosions. We plot in Fig. 
El the bolometric light curves for Models 2D-la-DDT us¬ 
ing the one-dimensional radiative transfer code. The nickel 
abundance and its distribution have stronger effects on the 
peak luminosity and width. We obtain a peak luminosity 
logio Tpeak = 43.1 with a drop of the absolute magnitude at 
15 days after maximum Amis = 0.58. A mild but observ¬ 
able secondary bump can be seen at about Day 30, show¬ 
ing that the photosphere has receded inside the iron layer. 
Beyond Day 40, when 56 Co decay becomes the dominant 
energy source, the curve drops at a constant rate. 


the bolometric light curve. We have preformed benchmark 
tests in one and two dimensions to evaluate the accuracy 
and robustness of the code. Also, we have studied the SNIa 
explosion using the PTD and DDT models. The nucleosyn¬ 
thesis details and the bolometric light curves of our sim¬ 
ulation models are comparable to those similar explosion 
models presented in the literature. 

However, it should be noted that a consistent compari¬ 
son with models in other works are difficult for two reasons. 
First, there are different numerical treatments for the same 
physical process. For example, in earlier works the effective 
flame propagation speed is chosen to be the maximum of 
various flame propagation mechanisms, while in later works 
the transition of the fla me fro m laminar to turbulent is based 
on analytic models (jPocheaul . ll994l l. However, how the flame 
propagation varies with local turbulence is not yet exactly 
known. Second, each code has its own prescription for the 
input physics, such as the EOS and the nuclear reaction net¬ 
work, while some processes, such as the nucleosynthesis and 
the turbulence model, depend on these m odels. For exam¬ 
ple, w e use the same EOS presented in Timme s fe Swestvl 
(Il999h as the FLASH code (ICalder et al.l . 2002fl . but their 
model does not include the effects of sub-grid turbulence. 
Also, our nuclear energy production depends on the choice 
of EOS. We directly use a nuclear reaction network to calcu¬ 
late the energy released by the deflag ration and detonation 
waves, while in ISchmidt et al.l (120051 1 all materials are as¬ 
sumed to be burnt into nuclear statistical equilibrium at 
high densities and only up to 24 Mg at low densities. The 
difference has a subtle effect on the explosion since the pro¬ 
duced energy may boost the fluid motion for the creation of 
sub-grid velocity fluctuations and at the same time makes 
the star expand. The produced energy then affects the local 
density which consequently affects the energy release again. 
Therefore it is difficult to carry out a quantitative compari¬ 
son when the hydrodynamics, flame physics and nucleosyn¬ 
thesis are coupled. 

In future work, we plan to include other numerical 
components, such as the advection-diffusion-reaction scheme 
and other models of sub-grid turbulence. By comparing the 
results of these numerical components we may better under¬ 
stand their roles in the explosion, as well as the limitation 
and validity regimes of these models . Also, as a nat ural con¬ 
tinuation of our previous work dLeung et all 12013I) . the ef¬ 
fects of dark matter will be included into the dynamics. We 
plan to examine whether the presence of dark matter affects 
the role of SNIa as a standard candle, which is critical in the 
discovery of dark energy and measurements of cosmological 
distances. 


7 DISCUSSION AND CONCLUSION 

In this paper, we have presented a two-dimensional Eule- 
rian hydrodynamics code for modeling SNIa using cylin¬ 
drical coordinates with the level-set method as the flame 
capturing scheme. We also included the sub-grid turbulence 
for modeling turbulent flame and a neutrino subroutine for 
computing the corresponding neutrino spectra. We also used 
a tracer particle scheme for nucleosynthesis and developed 
a moment-integrated radiative transfer code for computing 
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Table 7. Masses of all isotopes at the end of simulations for different numbers of tracer particles N = 40 2 , 80 2 and 160 2 . The approximate 
yields from the hydrodyanmics run are also included for comparison. The 19-isotope network is used. The record time is fixed at t = 0.025 
s. All quantities are in units of solar mass. 


Isotope 

Hydro 

N = 

40 2 

N = 

<N 

o 

00 

N = 160 2 

12 c 

6.68 x lO" 2 

1.61 x 

lO" 2 

1.61 x 

10 -2 

1.60 x 10 -2 

16 o 

6.68 x lCT 2 

3.74 x 

10“ 2 

4.36 x 

10“ 2 

4.33 x 10“ 2 

24 Mg 

1.60 x 10 -5 

8.36 x 

10 -3 

1.11 X 

10 -2 

1.06 x 10“ 2 

28 Si 

9.17 x 10" 3 

1.97 x 

10" 1 

1.99 x 

10" 1 

1.98 x 10" 1 

32g 

9.93 x 10" 3 

9.22 x 

10 -2 

9.16 x 

10 -2 

9.08 x 10" 2 

36 Ar 

6.90 x 10 -3 

1.93 x 

10 -2 

1.90 x 

10 -2 

1.88 x 10“ 2 

40 Ca 

1.27 x 10" 2 

1.82 x 

10 -2 

1.78 x 

10 -2 

1.76 x 10" 2 

44 Ti 

3.80 x 10" 4 

9.96 x 

10" 5 

9.73 x 

10" 5 

9.66 x 10" 5 

48 Cr 

2.91 x 10~ 2 

8.71 x 

10- 4 

8.46 x 

10- 4 

8.39 x lO" 4 

52 Fe 

4.20 x 10" 2 

2.10 x 

10“ 2 

2.07 x 

10 -2 

2.06 x 10 -2 

54 Fe 

0.368 

0.100 

9.89 x 

10 -2 

9.99 x 10" 2 

56 Ni 

0.733 

0.836 

0.831 

0.831 


APPENDIX A: DEFLAGRATION WAVE 

In Sections m and EJ we mentioned that deflagration and 
detonation are coupled to the hydrodynamics as energy 
sources. Here we describe the governing equations for de¬ 
flagration, how they are solved to obtain the results. 

Deflagration is one of the combustion processes where 
electron scattering is the major heat transport mechanism. 
The process is so slow that the fluid remains almost iso- 
ba ric. The wave struct u re an d the flame speed are studied 
in iTimmes fe Wooslevl dl992l l. At least four methods can 
be used to study the structure of deflagration wave. They 
include direct simulations, diffusion approximations, eigen¬ 
value method and variational approximations, arranged in 
descending order of information provided. Direct simulation 
is the computationally most expensive method, but is the 
most accurate one, because no assumption is made in the 
evolution equations. On the other hand, variational approx¬ 
imation is the least accurate as it gives only a lower bound of 
the flame speed and reveals no information about the flame 
structure. In our simulations, we use the diffusion approxi¬ 
mations in finding the flame structure and its corresponding 
speed, which is shown to be a good approximation due to 
the slow propagation of flame compared with the much faster 
sound speed. 

The diffusion approximation starts from the general Eu¬ 
ler equations in spherical coordinate and Lagrangian formu¬ 
lation, namely, 


dm 

dr 

= Anr 2 p, 

(Al) 

Dv 

~Dt 

„ , dP GM 

= —Anr - - —, 

am r z 

(A2) 

De d(l/p) 

Dt dt 

1 d ( dT\ . 

(A3) 

a 

- 

(A4) 

DYi 

Dt 

= I] -YiYjXjkii) + YAkXkjil). 

(A5) 


j,k 


Bi and Yi are the binding energy and the number fraction 
of isotope i. Xjk is the reaction rate of nuclear reaction from 


isotope j to k. Other variables have the the same physical 
meaning as in the main text. The derivative D/Dt is the 
time derivative in the frame moving with the fluid. The dif¬ 
fusion approximation then assumes that the fluid has always 
the same pressure. Gravity is neglected because the typical 
size of deflagration reaction zone is much smaller than the 
density scale height. This means that the velocity equation 
Eq. (IA2I) can be neglected. 

In the heat diffusion term on the right hand side of Eq. 
<lA3l) X is the total thermal conductivity, which consists of 
contributions from both electron and photon conductivities 

X = Xe + X 7 - ( A6 ) 

Photon conductivity can be exactly found as x-y = 
AacT 3 /3Kp, with a = 7.5657 x 1CP 15 erg cm -1 K -4 be¬ 
ing the radiation constant, c being the speed of light and re 
being the opacity. 

To find the electron condu ctivity, we follow the prescrip¬ 
tion of lKhokhlov et all (Il997l'l . First we have \ e expressed 
in terms of the effective electron collisional frequency 
3 /in 16 A 

Xe = 4.99x1 n 9 T — ergs cm -1 s -1 K -1 , (A7) 

S/l + X 2 \ Ve J 

with x = pF/m-eC being the dimensionless Fermi momen¬ 
tum. The electron collisional frequency can be separated into 
ion-electron and electron-electron collisional frequencies by 

V e =IS ee + Vei- (A8) 


The electron-ion collisional frequency is given by 

_ 7 A 

u Bi = 1.78 x 10 16 ^/l + x 2 ^—, (A9) 

AY e 

with A being the mean atomic mass and Z being the mean 
atomic number. Y e is the electron fraction and A is the 
Coulomb logarithm. 


A 


In 


^y /3 ^/ 3r + 3 


2(1 +a; 2 ) + 


—a/3 2 _ 1 + L3q _ 

2 M 1 + a 2 (0.71 — 0.54/3 2 ) 

Here, F = 2.275 x 10 5 Z 5 / 3 (pT e ) 1 / 3 /T, /3 
a = Z/137/3. 


(A10) 

= x/\/l + x 2 and 
























18 S.-C. Leung et al. 



Figure Al. Isotope abundances against density for 13 isotopes. 
See the text for the lists of elements. 



Figure A2. Isotope abundances against density for 19 isotopes. 
See the text for the lists of elements. 

Similarly, we have the electron-electron collision fre¬ 
quency as a function of temperature and density 

. x 3/2 

^ = °- 511T (1 + 22 ) 5/4 J (») 8 > ( A11 ) 

with y = \/3T pe /T. T pe is the electron-plasma temperature 
given by 

x 3/2 

Tpe = 3-307 x 10 8 (1 + x2)1/4 K. (A12) 

J(y) is an integral that cannot be evaluated analytically and 
is given numerically by 



where a = 0.113 and b = 1.247. 

Using these information, we find numerically the defla¬ 
gration wave propagation speed, the energy production and 
the ash composition. T o start the d eflagrati on wa ve, we fol¬ 
low the prescription in lTimmes fc Wooslevl (j l992h . We first 
ignite the innermost grid cells, setting its initial temperature 
sufficient for nuclear runaway. Then we start the evolution 



Figure A3. Specific internal energy production against density 
for 13, 19 and 489 isotopes. See the text for the lists of elements. 

and gradually a steady deflagration wave forms. The lami¬ 
nar flame velocity can be obtained by studying the motion 
of the steady wave structure. The ash composition is ob¬ 
tained from the regions swept by the deflagration wave. The 
internal energy production is calculated by comparing the 
initial and final sp ecific internal energy . _ 

As shown in iTimmes fc Wooslevl (Il992h . the laminar 
flame velocity can be well approximated by 

ui am = 92.0kms(of) , (A14) 

with Xc being the mass fraction of 12 C. 

We plot in Figs. lAll and lA2l the ash composition against 
fuel density using 13-isotope and 19-isotope networks. At 
low density p ~ 10 7 g cm -3 , only carbon is consumed. Oxy¬ 
gen remains barely unchanged. Magnesium and silicon are 
the major products. At p ~ 5 x 10' g cm -3 , oxygen is 
also consumed, with silicon and sulphur being the major 
products. At density around 10 s g cm -3 , iron and nickel 
become prominent. Also helium is produced due to photo¬ 
disintegration. 54 Fe is the major product when it is included 
in the nuclear reaction network, which is absent in the 13- 
isotope network. 

In Fig. lA3l we plot the energy release against fuel density 
using 13-, 19- and 489-isotope networks. The specific internal 
energy production (SIEP) is in the order of 10 17 erg. At low 
density p ~ 10 7 g cm~ 3 , SIEP increases with density and the 
oxygen mass fraction decreases. This is because the reaction 
rate increases with density. At density above 5 x 10' g cm' 3 , 
SIEP drops. This does not mean les s fuel is consumed. As 
mentioned in iReinecke et all (l2002al h at high temperature, 
the highly endothermic photo-disintegration 56 Ni —» 14 4 He 
becomes important, which compensates for the rise of SIEP. 


APPENDIX B: DETONATION WAVE 

To d eterm ine the wave structure, we follow the prescription 
in HE arpel (|l999h . We first find the thermodynamics state of 
shocked fluid from a given state. Then, using the shocked 
state as an initial condition, the wave structure can be ob¬ 
tained by solving the steady state limit of Euler equations 
with nuclear reactions. 
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For the first part, the upstream and downstream states 
are related by the Rankine-Hugoniot conditions 


pu = p 0 D, 


i 2 

P + pv 



= Po + poD 2 , 


= eo + 


Po 

po 



(Bl) 

(B2) 

(B3) 


Quantities with a subscript 0 are the upstream (unshocked) 
state. D is the outflow velocity. Both eo(po,To,Xo) and 
p(po,To, Xo) are dependent on given density, temperature 
and composition. In particular, po and To are free parame¬ 
ters, but the actual results are insensitive to To, because the 
upstream electrons are highly degenerate. Xo is the compo¬ 
sition of fuel, namely 0.5 X'c and 0.5 Xo- In general, the 
final composition X is a set of quantities (depending on the 
choice of isotopes) to be determined. 

After having the post-shock state, we determine the re¬ 
action zone size and the ash composition by assuming a one¬ 
dimensional planar and steady flow, namely 


dp du 

V te +P d^ = °’ 

(B4) 

du dp 

pv Tx + Tx=^ 

(B5) 

d± _ P_dp _ 

dx p 2 dx ’ 

(B6) 

dX R 

(B7) 

dx v ’ 


where R is the reaction rate. Then, using 


*-(&)++(£) .."■+£ 


dp 


dp 


(dp_ 


dT Kx V dXi J 

/ P1=1 x 7 Pi 


dXi , 
(B8) 


and 
dt = 


dt\ 

dp) 


) dp + 
T.X 


dt_ 

dT 


p,x 


dT +J2 


(— 

\dXi 


P,T 


dXi, 
(B9) 

the Euler equations with nuclear reactions in the steady 
state limit can be reduced to 


dp 

dx 

dT 

dx 


dY 

dx 

where 


pCLf a ■ R 

V 7] 

dpy 1 
d T ) p , x 

N 

E 


T,X 


dp 

dx 


dp 

~dXi 


P.T,X jfti 


dXj 

dx 


R 

v ’ 


2 2 
7] = Of — V 


(BIO) 

(BH) 

(B12) 

(B13) 


is the sonic parameter, 




dp 


dp): 


dp 

dT 


p,x 


dt_ 

dT 


p, x 
(B14) 

is the sound speed of constant composition (also known as 


frozen sound speed) and 


Oj 


1 


dp 

dXi 


P.T.X&i 





(B15) 


is the thermicity constant, such that a • R is the thermicity. 

Eq. (IB12I) can be integrated into the reaction zone. No¬ 
tice that there are actually two types of detonation struc¬ 
ture implied from these equations. One type is that through¬ 
out the detonation wave both p and thermicity are positive, 
while the other type is that both quantities are both positive 
or negative simultaneously. The first type is the well-known 
Chapman-Jouguet (CJ) detonation, while the second type 
is called the supported pathological (SP) detonation. 

CJ detonation occurs at low density (p < 2 x 10 7 g 
cm' 3 ). The boundary condition is given by p — a ■ R = 0 at 
x —> oo. This corresponds to the ash propagating at frozen 
sound speed at the point where no more nuclear reaction 
can proceed. 

For SP detonation (p > 2 x 10' g cm -3 ), the integration 
is divided into two sections. First, we integrate the equation 
close to p —> 0 (while thermicity needs not to be zero at that 
position). Assume p = 0 at x = x p and we have integrated 
up to i = - Xx p , by making use of the symmetry of 

Eqs. (1B12I) from the zero-velocity point and the continuity 
of thermodynamics variables, we obtain the post-zero-point 
states at x — x p + Ax p , where 


p(x p + Ax p ) = p(x p — Ax p ), (B16) 

and 


T(x p + Axp) = T(x p — Axp) - 

N 


E( — 

^ \ dXi 


P.T.X&, 


^r(2Ax p ). (B17) 


After the above transition, the integration can be carried on 
again until no net nuclear reaction continues. To determine 
the correct eigenvalue, we require p = a ■ R = 0 at the same 
position. 


APPENDIX C: FLAME CAPTURING USING 
THE POINT-SET METHOD 


In two-dimensional simulations, the flame front is repre- 


(see for example in (Glimm 

et al., 

1981 

Glimm & McBrvan, 

19851; Glimm et all 1987, 

19881, 

20 (L 

) for applications in 

two-dimensional systems, 

1 Glimm et all Il99fll 200C 

, 2003f) 


f or a pplications in three-dimensional systems and Zhanei . 
l2009h for its application in SNIa modeling.) introduces a set 
of pseudo-particles, which form a line by assuming that the 
nearest neighbors are linked. Similar to the level-set method, 
the particles are transported by the fluid advection with a 
speed Ffluid and its own propagation with velocity Uflam B n , 
such that 


Unode i — ^flam &7li T Ufl u id- (Cl) 

hi is the unit normal vector of the flame front pointing from 
the i-th node towards the fuel. The fluid velocity is given 
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directly by the Euler equations. The flame propagation de¬ 
pends on the fluid density and the local normal direction of 
flame front pointing towards the fuel. The normal of a node 
is defined by the positions of its closest neighboring nodes. 
For example, given a node i with position with a 

distance to its next node i + 1 given by 

j-1 — (Xj Xi-^. l)^ T ( yi Vi-\- 1)^; (C2) 


we define the normal direction of the line segment formed 
by the nodes i and i + 1 by 


f 2/i+l Vi *£i+1 
\ 1 1 


(C3) 


Notice that whether the normal vector pointing towards the 
fuel or ash is arbitrary. In our study, we choose the former 
one. 

After all node velocities are found, the positions are 
updated by 


■^node (new ) — 'T'node (old) T Unode A£ no de • ( C 1) 

There is no accelerating term for the node because the parti¬ 
cles serve only as markers of the deflagration front for com¬ 
puting the change in energy and isotopes, which do not in¬ 
teract with the fluid directly. 

The time step may not be the same as that of the hy¬ 
drodynamics one, because the resolution of the point-set 
method and the hydrodynamics are independent. In gen¬ 
eral, we require 

Af node = C-^-, (C5) 

^node 

with C a positive number smaller than unity and Z m i n the 
minimum distance between neighboring nodes. 

Apart from Z m i n , there are two more parameters control¬ 
ling the point-set resolution, Z max and Z me rge- They define 
the inter-node maximum separation and merging distance 
between non-neighboring nodes. These parameters are es¬ 
sential in maintaining a consistent resolution of the surface 
during the evolution, because the distance between linked 
nodes can become too far or close, or there can be lines cross¬ 
ing each other. These phenomena occur frequently when 
the fluid motion is turbulent. In those cases, node addition, 
removal or reconnection is needed. Notice that in the lit¬ 
erature, node reconnection is also known as surface split¬ 
ting/merging in the context of multi-dimensional simula¬ 
tions. 

To check whether there are regions which are over¬ 
crowded or underpopulated with nodes, we calculate the 
distance between nearest neighbors. Given two nodes i and 
i +1 with separation cZ;,; + i, if di^+i < Z m i n , one of the node is 
moved to a new position ((Xi + x i+ 1 )/ 2 , {yi + Vi+ 1 )/ 2 ), while 
the other node is deleted; on the other hand, if di ,i+1 > l max 
an extra node between node i and i + 1 is inserted at posi¬ 
tion ((*i-|- 2 ;i + i)/ 2 , (yi+y i+ i)/2). See Fig.[CT]for a graphical 
illustration of the above operations. 

To avoid lines from crossing each other, we locate all 
node pairs which are not connected but are potentially close 
enough to form a line. We first identify the nodes which lie 
on the same or neighboring grids, based on the same Eu- 
lerian grid of the hydrodynamics, then the separations of 
this group of nodes are computed. When any pair of non¬ 
neighboring nodes satisfies di,j < Z merge for some i and j that 
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Figure Cl. Illustration on node addition, removal and reconnec¬ 
tion. The arrows between nodes stand for the pointer from one 
node variable to the next node variable. 


|i — j | ^ 1, we reconnect the nodes as shown in Fig. IClI and 
change the topology of the surface. In principle, the aim is 
that the new surface may preserve the node position but the 
new surfaces will propagate away from each other, so that 
no entanglement can be formed in the coming time steps. In 
practice, for a line without entanglement, there exists only 
one unpaired node j with any node i such that dij < Z max . 
However, when surfaces are going to split or merge, there is 
more than one of such candidate, which means that there 
is more than one way to form a line. To decide which node 
should be chosen, we use the above principle to connect the 
node such that a concave surface is obtained. Geometrically, 
the normal vectors on the new surfaces are converging. See 
Fig. lC2l for a graphical description. We remark that how the 
surfaces merge is unique in two-dimensional simulations be¬ 
cause all nodes have at most two neighbors for forming line 
segments. This property is not true for three-dimensional 
models because the surface is usually represented by tri¬ 
angular patches, which means that each node always has 
multiple connections with neighboring nodes. Thus, the ge¬ 
ometry of the surface, such as the curvature, will depend on 
how nodes are connected. 

The three parameters, Z max , Z m i n and Z mer ge are inter¬ 
related. First, we require Z max = 2Z m i n . This is because we 
do not want any repetitive loop of node addition or removal 
to take place, that means, after adding (removing) a new 
node between any two nearest-neighbors, the new configu¬ 
ration does not have nodes which are too close to (far from) 
each other according to the criteria. Similarly, we require 
Zmerge < Z max to avoid entangled lines in simulations. Since 
the point-set method has a time step that prevents a node 
from moving further than Z max , we thus choose Z merg e = Z max . 

We compare the performance of the level-set algorithm 
with the point-set algorithm. The level-set algorithm is 
known to be unable to track the flame consistently when the 
flame propagates much slower than the fluid flow. The flame 
surface is dominated by the turbulent flow where the flame 
is supposed to show a convoluted and elongated structure 
shaped by the fluid. But the level-set method can only pre¬ 
serve structure with a size greater than the grid size. There¬ 
fore, the small-scale structure cannot be tracked and a de- 
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After merging 


Figure C2. Illustration on how node merging can pervent two 
lines from crossing each other. The arrows on the line stand for 
the normal direction of the flame surface. 



r (km) 

Figure C4. Flame front of Model PSM-1-2 at t = 1.00 s. 



time (s) 


Figure C3. Total energy against time for Models LSM-1, PSM- 
1-2, PSM-1-3, PSM-1-4. See Table ICll for the configurations. 



r (km) 


Figure C5. Flame front of Model PSM-1-3 at t = 1.00 s. 


tached flame in forms of bubbles is obtained. To maintain the 
consistency of the level-set algorithm in SNI a simu l ations , 
flame acceleration schemes are needed (ICalder et all [200711 . 
To address this problem, we compare the flame structure 
in the laminar flame limit by using the two algorithms. We 
present the simulation models in Table E2 The hydrodyan- 
mics is the same as those in PTD tests and DDT tests. The 
hydrodynamics is done with a configuration similar to Sec¬ 
tion [5]3] with an array of 500 x 500 in cylindrical coordinates 
with uniform grid size A = 11 km. 

In Fig. IC3I we plot the total energy against time for 
Models LSM-1, PSM-1-2, PSM-1-3 and PSM-1-4. As ex¬ 
pected, the laminar deflagration cannot successfully unbind 
the star. All four models perform similarly in early time. The 
model PSM-1-3 is the most similar one to LSM-1. At later 
time, the level-set method predicts less energy as compared 
to other three models. 

In Figs. IC4IIC5IIC6I and IC7I we plot the flame front of 
Models PSM-1-2, PSM-1-3, PSM-1-4 and LSM-1, respec¬ 
tively at t = 1.0 s. By comparing Models PSM-1-2, PSM- 
1-3 and PSM-1-4, it shows that the resolution of point-set 
method affects the performance in two ways: First, it con¬ 


tributes to the fine details of RT and KH instabilities, which 
enlarge the local flame area; second, the algorithm allows an 
extremely elongated flame structure and an injection of fuel 
into burnt region in the sub-grid scale. This property is im¬ 
portant for keeping information of the flame surface as much 
as possible, which maintains the fuel burning rate. 

We compare the point-set method and the level-set 
method by comparing the flame surface of LSM-1 and PSM- 
1-4. Both methods give thin flame shape. The point-set 
method shows a largely connected structure with some fuel 
regions surrounded by ash. In level-set method, the flame is 
broken into pieces and the flame bubbles are disassociated 
and are apart from each other. There is no sign of fuel in¬ 
jection. The upper part of the flame can be compared with 
Model PSM-1-2. But the narrow parts in Model PSM-1-2 
which extends from the core are not seen in LSM-1. Only 
the largest structure above the resolution size is preserved. 
Comparing these results, the point-set method can success¬ 
fully capture the flame structure even in the laminar flame 
limit, which is essential in understanding the energy release 
rate of highly convoluted flame. 
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Table Cl. Simulation setup for the test of flame capturing scheme. Lengths are in unit of code unit and densities are in unit of 10® g 
cm -3 . 


Model 

scheme 

/max 

/min 

/merge 

Pc 

LSM-1 

level-set 

3.0 

1.5 

3.0 

3 

PSM-1-2 

point-set 

2.0 

1.0 

2.0 

3 

PSM-1-3 

point-set 

3.0 

1.5 

3.0 

3 

PSM-1-4 

point-set 

4.0 

2.0 

4.0 

3 



r (km) 

Figure C6. Flame front of Model PSM-1-4 at t = 1.00 s. 



r (code unit) 


Figure C7. Flame front of Model LSM-1 at t = 1.00 s. 


APPENDIX D: IONIZATION AND OPACITIES 


We use an open-source Saha-equation solver0. The sub¬ 
routine solves the Saha equations, which relate the number 
densities of the i th -times ionized and the i + 1-times ionized 
species by 


ni+i t zn e 

n it z 




(Dl) 


where t&gz is a function of the partition functions of both 
species at a given temperature and ionization energies. In 
the subroutine, elements from 1 H to 30 Zn are included and 
all ionization stages are considered. An initial guess is given 
and then the solution is obtained by iterations. 

The ionization fractions of all elements are then 
summed to find the number density of free electrons, which 
is used by the Helmholtz EOS subroutine in order to solve 
for the hydrodynamics pressure, internal energy and other 
local thermodynamics quantities. 

After the number densities of free electron n e and the 
number fraction nz,i of an element Z in an ionization stage i 
are found, we obtain the total opacity k, which includes the 
Thomson opacity K e , bound-free opacity K b ? (u) and free-free 
opacity K bb (p) 

k{v) = K e + K b f (v) + K bb (v). (D2) 


The Thomson opacity is given by 


K, e 


TLeCT e 

P 


(D3) 


where cr e ~ 6.65 x 10~ 25 cm 2 is the Thomson cross-section of 
electron. The Thomson opacity is temperature and density 
independent to a good approximation. 

The bound-free opacity, or the photo-ionization opacity, 
is given by 


bf / \ a z,i( l/ ) nz ’ i 

« H = v —- 

' p 

Z,i y 


(D4) 


with cr^)(i/) being the frequency-dependent bound-free scat¬ 
tering cross-section for the element Z at the i th ionization 
stage. We use the fitting formula reported bv IVerner et al.l 
d 19961) 


In Section l4~2l we described a hydrodynamics scheme in mod¬ 
eling the evolution of radiation of matter for a few weeks 
after the SNIa explosion. In general, due to the rapid ex¬ 
pansion and radiation, the matter reaches a temperature 
where the assumption of complete ionization becomes in¬ 
valid. Therefore, including ionization fraction in the light 
curve modeling is important for a consistent description of 
the thermodynamics properties of the matter. 


az,i( v ) = <ro - i ) 2 + yw\ y p/2 5 ' 5 ( 1 + \/^f) P > 

V i/a 

(D5) 

2 Provided as open - source code in www.cococubed.edu. Refer 
llPaxton et all 1201(1 l2013h for its applications in stellar evolu¬ 
tion. 
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where 

x{v) = — - yo (D6) 

eo 

and 

V = \Jx 2 + yf (D7) 

The constants eo, 00 , yo, yi, ya, y-w and P are constants 

depending on the species and ionizat ion stage. _ 

The free-free opacity is given by dSakamoto et alJ . 1200Ch 

k !S { v) = Cz 1 in e nz,iP'/Ti'- S ', (D8) 

Z,i 


with Cz,i being some const ants related to the Gaunt factor , 
whose values are fitted in (lltoh et all 1 19861 : iNozawa et all 
1 19961: lltoh et all l2000ll . 

The Rosseland mean opacity is then given by 


Rr 


r 00 1 dB u (T) 

Jo K(i 7 ) dT 


du 


r 00 dB v (T) 


JO 


dT 


dv 


(D9) 


with B„(T) being the Planck distribution function. Since 
only the bolometric light curve is modeled, we assume Kp = 
ke and = XR i n the calculation. 

The effects of lines are important because of the Doppler 
effect inside the fast expanding ejecta. Photons with fre¬ 
quencies within Doppler widths from the atomic scatter¬ 
ing lines can be absorbed or scatter ed. To include the line 
opacity, we fol l ow th e prescription in lKarp et al.l dl977l 'l and 
Hocfl ich et al.l (Il993l h We first compute the Sobolev opacity 
of each line by 


^line(f, j) 


ne 2 fijm / _ gjUj 
rrieC Uijp \ gjm 


(DIO) 


with fij being the oscillator strength between atomic states 
i and j, rii and gi are the occupation number and statistical 
weight at state i. Vij is the transition frequency. 

The effective absorption coefficient is obtained by 
summing all relevant lines, where the lines are arranged in 
increasing wavelength 


X ^ — ^cont X [1 

Y 1 - exp(-Tj) ( } exp 

j=J V 

with 

dr 

Ti — KUneCp— r 

dv(r) 



(D12) 


being the effective optical depth of the line i. Kcont is the 
continuum opacity (bound-free opacity, free-free opacity and 
Thomson scattering). The frequency dependent expansion 
rate parameter s{y) is given by 

s{y) = Kcont cp (D13) 

We remark that as pointed out in iBlinnikovI d 1996h and 
IPinto et al.l d2000l l. the summation in Eq. ID11I should be 
done in an upwind direction so that the opacity includes 
only the lines that are being scattered. 

To relate the extinction coefficient with the absorption 


coefficient, we define the frequency independent enhance¬ 
ment factor 


e 



11 - 


with xr the same Rosseland opacity but with the line opac¬ 
ity included. To discriminate the difference between line ab¬ 
sorption and line scattering, we define the ratio between line 
collisional and radiative transition rate between state i and 
j by q, which modifies the the enhancement factor by an 
extra factor 


a 

a + l 6 ' 


(D15) 


The frequency dependent absorption coefficient is 


Kv = e{Xv - CTe). (D16) 


References 

Arnett W. D., 1969, Ap & SS, 5, 180 
Arnett W. D., 1982, Astrophys. J., 253, 785 
Barth T. J., Deconinck H., 1999, Lecture Notes in Compu¬ 
tational Science and Engineering 9: High-Order Methods 
for Computational Physics, Springer. 

Bell J. B., Day M. S., Rendleman C. A., Woosley S. E., 
Zingale M., 2004, Astrophys. J., 606, 1029 
Blinnikov S. I., 1996, AZh Pisma, 22, 92 
Blinnikov S. I., Sorokina E. I., 2000, A & A, 356, L30 
Blinnikov S. I., Eastman R., Bartunov O. S., Popolitov 
V. A., Woosley S. E., 2008, Astrophys. J., 496, 454 
Blinnikov S. I. et al., 2006, Astrophys. J., 453, 229 
Blondin S., Dessart L., Hillier D. J., Khokhlov A. M., 2012, 
MNRAS, 429, 2127. 

Blondin S., Kasan D., Roepke F. K., Kirshner R. P., Man- 
del K. S., 2011, MNRAS, 417, 1280 
Branch D., Tammann T. A., 1992, Ann. Rev. Astron. As¬ 
trophys., 30, 359 

Calder A. C. et al., 2002, Astrophys. J. Suppl., 143, 201 
Calder A. C. et al., 2007, Astrophys. J., 656, 313 
Clement M. J., 1993, Astrophys. J., 406, 651 
Colgate S. A., Petschek A. G., Kriese J. T., 1980, Astro¬ 
phys. J., 237, L81 

Colgate S. A., White R. H., 1966, Astrophys. J., 143, 626 
Collela P., Woodward P. R., 1984, J. Comput. Phys., 54, 
174 

Damkoehler G., 1939, Jahrb. Deut. Luftfahrtforsch., 113 
Dessart L., Blondin S., Hillier D. J., Khokhlov A. M., 2013, 
441, 542. 

Timmes F. X., Arnett D., 1999, Astrophys. J., 125, 277 
Ferreira V. G., et al., 2004, Int. J. Numer. Meth. Fluids, 
44, 347 

Fink M., et al., 2014, Mon. Not. R. astr. Soc., 438, 438 
Fryxell et al., 2000, ApJS, 131, 273 

Gamezo V. N., Khokhlov A. M., Oran E. S., 2004, Phys. 
Rev. Lett., 92, 211102 

Gamezo V. N., Khokhlov A. M., Oran E. S., 2005, Astro¬ 
phys. J., 623, 337 

Glimm J., Isaacson E., Marchesin D., McBryan O., 1981, 
Advances in Applied Mathematics, 2, 91 
Glimm J., McBryan O. A., 1985, Advances in Applied 
Mathematics, 6, 422 










































24 S.-C. Leung et al. 


Glimm J., McBryan O., Menikoff R., Sharp D. H., 1987, 
SIAM J. Sci. Stat. Comput., 7, 230 
Glimm J., Grove J., Lindquist B., McBryan O. A., Tryg- 
gvason G., 1988, SIAM J. Sci. Stat. Comput., 9, 006 
Glimm J., Grove J. W., Li X. C., Zhao N., 1999, Contem. 
Maths., 238, 133 

Glimm J., Grove J. W., Li X. C., Tan D. C., 2000, SIAM 
J. Sci. Comput., 21, 2240 

Glimm J., Grove J. W., Zhang Y., 2002, SIAM J. Sci. Com¬ 
put., 24, 208 

Glimm J., Li X., Liu Y., Xu Z., Zhao N., 2003, SIAM J. 
Num. Anal., 41, 1926 

Golombek I., Niemeyer J. C., Astron. Astrophys., 2005, 
438, 611 

Gresho P., Int. J. for Num. Meth. in Fluids, 1990, 11, 621 
Gueyffier D., Li J., Nadim A., Scardovelli R., Zaleski S., 
1998, J. Comput. Phys, 152, 423 
Hauschildt P. H., Baron E., 2010, A & A, 509, A36 
Hauschildt P. H., Baron E., 2011, A & A, 533, A127 
Hillier D. J., 1990, A & A, 231, 111 
Hillier D. J., 1990, A & A, 231, 116 
Hillier D. J., Miller D. L., 1998, ApJ, 496, 407 
Hirt C. W., Nichols B. D., 1981, J. Comput. Phys., 39, 201 
Hoeflich P., Mueller E., Khokhlov A., 1993, Astron. Astro¬ 
phys., 268, 570 

Hoeflich P., Khokhlov A., Wheeler J. C., 1995, ApJ, 444, 
831 

Imshenik V. S., Kal’yanova N. L., Koldoba A. V., 
Chechetkin V. M., 1999, Astronomy Letters, 25, 206 
Itoh N., Nakagawa M., Kohyama Y., 1986, ApJ, 294, 17 
Itoh N., Hayashi H., Nishikawa A., Kohyama Y., 1996, ApJ, 
102, 411 

Itoh N., Sakamoto K., Kusano S., Nozawa S., Kohyama Y., 
2000, ApJS, 128, 125 

Jiang G.-S., Shu C.-W., 1996, J. Comput. Phys., 126, 202 
Jorden IV G. C. et al., 2008, ApJ, 681, 1448 
Jorden IV G. C. et al., 2012, ApJ, 759, 53 
Karp A. H., Lasher G., Chan K. L., Salpeter E. E., 1977, 
ApJ, 214, 161 

Kasen D., 2006, Astrophys. J., 651, 366 
Kasen D., Plewa T., 2007, Astrophys. J., 662, 459 
Khokhlov A. M., 1989, Mon. Not. R. astr. Soc., 239, 785 
Khokhlov A., 1993, Astrophys. J., 419, L77 
Khokhlov A., 1994, Astrophys. J., 424, L115 
Khokhlov A. M., 1991a, Astron. Astrophys., 245, 114 
Khokhlov A. M., 1991b, Astron. Astrophys., 245, L25 
Khokhlov A. M., 1991c, Astrophys. J., 246, 383 
Khokhlov A. M., Oran E. M., Wheeler J. C., 1995, Astro¬ 
phys. J., 478, 678 

Khokhlov A. M., Oran E. S., Wheeler J. C., 1997, Astro¬ 
phys. J., 478, 678 

Kromer M., Sim S. A., 2009, Mon. Not. R. astr. Soc., 398, 
1809 

Launder B. E., Spalding D. B., 1974, Int. J. Numer. Meth. 
Fluids, 15, 127 

Leibundgut B., Pinto P. A., 1992, Astrophys. J., 401, 49 
Leung S.-C., Chu M.-C., Lin L.-M., Wong K.-W., 2013, 
Phys. Rev. D, 87, 123506 

Leung S.-C., Chu M.-C., Lin L.-M., in preparations. 
Lisewski A. M., Hillebrandt W., Woosley S. E., 2000, As¬ 
trophys. J., 538, 831 


Livne E., Ashida S. M., Hoeflich P., 2005, Astrophys. J., 
632, 443 

Long M. et al., 2014, ApJ, 789, 103 
Lucy L. B., 2001, MNRAS, 326, 95 
Lucy L. B., 2005, Astron. Astrophys., 429, 19 
Maier A., Niemeyer J. C., 2006, Astron. Astrophys., 451, 
207 

McNally C. P., Lyra W., Passy J. C., 2011, ApJS, 201, 18 
Meakin C. A. et al., 2009, ApJ, 693, 1188 
Misiaszek M., Odrzywolek A., Kutschera M., 2006, Phys. 
Rev. D, 74, 043006 

Mueller E., Arnett W. D., 1986, Astrophys. J., 307, 619 
Niemeyer J. C., 1999, Astrophys. J., 523, L57 
Niemeyer J. C., Bushe W. K., Ruetsch G. R., 1999, Astro¬ 
phys. J., 524, 290 

Niemeyer J. C., Hillebrandt W., 1995a, Astrophys. J., 452, 
769 

Niemeyer J. C., Hillebrandt W., 1995b, Astrophys. J., 452, 
779 

Niemeyer J. C., Hillebrandt W., Woosley S. E., 1996, As¬ 
trophys. J., 471, 903 

Niemeyer J. C., Woosley S. E., 1997, Astrophys. J., 475, 
740 

Nomoto K., Sugimoto D., 1977, Publ. Asron. Soc. Japan, 
29, 765 

Nomoto K., Sugimoto D., Neo S., 1976, Ap & SS, 39, 37 
Nomoto K., Thielemann F.-K., Yokoi K., 1984, Astrophys. 
J., 286, 644 

Nomoto K., Thielemann F.-K., Yokoi K., Branch D., 1986, 
Ap & SS, 118, 305 

Nozawa S., Itoh N., Kohyama Y., 1998, ApJ, 507, 630 
Nugent P., Baron E., Branch D., Fisher A., Hauschildt P. 
H., 1997, Astrophys. J., 485, 812 
Odrzywolek A., 2007, Eur. Phys. J. C, 52, 425 
Osher S., Sethian J. A., 1988, J. Comput. Phys., 79, 12 
Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., 
Timmes F., 2010, ApJS, 192, 3 
Paxton B. et al., 2013, ApJS, 208, 4 
Perlmutter S., et al., 1999, Astrophys. J., 517, 565 
Phillips M. M., et al., 1987, PASP, 99, 592 
Pinto P. A., Eastman R. G., 2000, ApJ, 530, 757 
Plewa T., 2007, Astrophys. J., 657, 942 
Pocheau A., 1994, Phys. Rev. E, 49, 1109 
Reinecke M., Hillebrandt W., Niemeyer J. C., 1999, Astron. 
Astrophys., 347, 739 

Reinecke M., Hillebrandt W., Niemeyer J. C., Klein R., 
Gloebl A., 1999, Astron. Astrophys., 347, 724 
Reinecke M., Niemeyer J. C., Hillebrandt W., 2002a, As¬ 
tron. Astrophys., 386, 936 

Reinecke M., Niemeyer J. C., Hillebrandt W., 2002b, As¬ 
tron. Astrophys., 391, 1167 

Rider W. J., Kothe D. B., 1995, AIAA Paper, 95, 1717 
Riess A. G., et al., 1998, Astron. J., 116, 1009 
Roepke F. K., 2005, Astron. Astrophys., 432, 969 
Roepke F. K., 2007, Astrophys. J., 668, 1103 
Roepke F. K., Hillebrandt W., 2005, Astron. Astrophys., 
431, 635 

Roepke F. K., Hillebrandt W., Niemeyer J. C., 2004a, As¬ 
tron. Astrophys., 421, 411 

Roepke F. K., Hillebrandt W., Niemeyer J. C., 2004b, As¬ 
tron. Astrophys., 421, 783 

Roepke F. K., Hillebrandt W., Schmidt W., Niemeyer J. C., 



A new hydrodynamics code for Type la Supernovae 25 


Blinnkov S. I., Mazzali P. A., 2007, Astrophys. J., 668, 
1132 

Roepke F. K., Hillebrandt W., Woosley S. E., 2007, Astro¬ 
phys. J., 660, 1344 

Roepke F. K., Niemeyer J. C., Hillebrandt W., 2003, As¬ 
trophys. J., 588, 952 

Rudman M., 1997, Int. J. Numer. Methods Fluids, 24, 671 
Sakamoto T., Itoh N., Kusano S., Nozawa S., Kohyama Y., 
2001, ASP Conf. Series, 251, 268 
Sauer D. N., Hoffmann T. L., Pauldrach, A. W. A., 2006, 
A & A, 229, 240 

Scadovelli R., Zaleski S., 1999, Ann. Rev. Fluid Mech., 31, 
567 

Schmidt W., Niemeyer J. C., Hillebrandt W., 2005, Astron. 
Astrophys., 450, 265 

Schmidt W., Niemeyer J. C., Hillebrndt W., Roepke F. K., 

2006, Astron. Astrophys., 450, 283 

Seelmann A. L., Hauschildt P. H., Baron E., 2010, A & A, 
522, A102 

Seitenzahl I. R., Roepke F. K., Fink M., Pakmor R., 2010, 
MNRAS, 407, 2297 

Seitenzahl I. R., et al., 2013, MNRAS, 429, 1156 
Sethian J. A., 1996, Proc. Natl. Acad. Sci. USA, 93, 1591 
Sethian J. A., 1999, SIAM Review, 41, 2 
Sethian J. A., 2001, J. Comput. Phys., 169, 503 
Sharpe G. J., 1999, MNRAS, 310, 1039 
Shih T.-H., Liou W. W., Shabbir A., Yang Z., Zhu J., 1994, 
Comp. Fluids, 24, 227 

Shih T.-H., Zhu J., Lumley J., 1995, Comp. Methods Appl. 

Mech. Engng., 125, 287 
Sim S. A., 2007, Astron. Astrophys., 375, 154 
Sim S. A., Sauer D. N., Roepke F. K., Hillebrandt W., 

2007, Astron. Astrophys., 378, 2 
Smagorinsky J., 1963, Mon. Weather Rev., 91, 99 
Strain J., 1999, J. Comput. Phys., 151, 498 

Sussman M., SmerekaP., Osher S., 1994, J. Comput. Phys., 
114, 146 

Swartz D. A., Sutherland P. G., Harkness R. P., 1995, As¬ 
trophys. J., 446, 766 

Timmes F. X., Woosley S. E., 1992, Astrophys. J., 396, 649 
Timmes F. X., 1992, Astrophys. J., 423, L131 
Timmes F. X., Hoffman R. D., Woosley S. E., 2000, Astro¬ 
phys. J. Suppl., 129, 377 

Timmes F. X., Swesty F. D., 1999, Astrophys. J. Suppl., 
126, 501 

Toro E., Riemann Solvers and Numerical Methods for Fluid 
Dynamics, Springer-Verlag, Berlin, Heidelberg, 1997. 
Towsley D. M. et al., 2007, Astrophys. J., 668, 1118 
Travaglio C., Hillebrant W., Reinecke M., Thielemann F.- 
K., 2004, A & A, 425, 1029 

Tryggvason G. et al., 2001, J. Comput. Phys., 169, 708 
Unverdi S. A., Tryggvason G., 1992, J. Comput. Phys., 100, 
25 

Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 
1996, Astrophys. J., 465, 487 

Wang R., Spiteri R. J., 2007, SIAM J. Numer. Anal., 45, 
1871 

Woosley S. E., Kasen D., S. Blinnikov, Sorokina E., 2009, 
Astrophys. J., 662, 487 

Woosley S. E., Kerstein A. R., Sankaran V., Aspden A. J., 
Roepke F. K., 2009, Astrophys. J., 704, 255 
Yoshizawa A., Abe H., Matsuo Y., Fujiwara H., Mizobuchi 


Y., 2012, Phys. Fluids, 24, 075109 
Youngren G. K., Acrivos A., 1976, J. Fluid Mech., 76, 433 
Zhang X., Sutherland P., 1994, Astrophys. J., 422, 719 
Zhang Y., 2009, Nonlinearity, 22, 1909 



